//------------------------------------------------------------------------------
//	AceFile.cc
//
//------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/AceFile.cc#3 $
//------------------------------------------------------------------------------
// $Log:$
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "AceFile.hh"

// I want reverse sort (descending), so define my own compare
bool Ace_lt_int(int c1,int c2){return(c2<c1);};


//---------------------------------------------------------------------------
// * FindToken
//---------------------------------------------------------------------------
// Returns true if tempLine contains the token string, scans through the file 
// looking for a line containing the token otherwise, returning true if found, 
// returns false if end of file is reached and the token isn't found. 
//
bool FindToken(istream &inFile,string feature,
	       string &tempLine)
{
  if (Contains(tempLine,feature)) return(true);  
  int idx;
  bool bFeatureMarker;
  do{
    GetLine(inFile,tempLine);
    idx = tempLine.find(feature);
    bFeatureMarker = (idx>=0);
  }while((!inFile.eof()) && (!bFeatureMarker));
  
  if (!bFeatureMarker){ 
    return(false);
  }else{
    return(true);
  }
}

//---------------------------------------------------------------------------
// * ContainsWSPadded
//---------------------------------------------------------------------------
// Utility function. 
// Returns true if the given feature appears in the query padded by white space. 
//

bool ContainsWSPadded(string &query,string feature)
{
  int idx = query.find(feature);
  if (idx >=0){
    if ((idx+feature.size()) == query.size()) return(true); // it's only thing on line. 

    int fsize = feature.size();
    if ((idx+fsize) < query.size()){
      if ((query[idx+fsize] == '\t') || 
	  (query[idx+fsize] == ' ') || 
	  (query[idx+fsize] =='\n')){
	return(true);
      }else return(false);
    }
  }else return(false);
}



//---------------------------------------------------------------------------
// * FindTokenPlusWS
//---------------------------------------------------------------------------
// Differs from FindToken in that it looks for Token followed by some white
// space.  This is to help avoid spurious matches when, for example, a tracename
// includes the Token  e.g.   TMBQ1U0001   is not to be confused with a line
// containing just BQ, or with a line containing BQ 2456  or some such.  
//
bool FindTokenPlusWS(istream &inFile,string feature,
		     string &tempLine)
{
  if (Contains(tempLine,feature)) return(true);
  
  int idx;
  bool bFeatureMarker;
  do{
    GetLine(inFile,tempLine);
    idx = tempLine.find(feature);
    bFeatureMarker = (idx>=0);
  }while((!inFile.eof()) && (!bFeatureMarker));

  int fsize = feature.size();

  if (bFeatureMarker){
    // check to see if this feature is it for this line, or if it is followed by white space. 
    if ((idx+feature.size()) == tempLine.size()) return(true);
    if ((idx+fsize) < tempLine.size()){
      if ((tempLine[idx+fsize] == '\t') || (tempLine[idx+fsize] == ' ') || 
	  (tempLine[idx+fsize] == '\n')) return(true);
    }
    cerr<<"Find Token How did we get here? "<<endl;
    cerr<<"feature="<<feature<<"idx="<<idx<<" linesize="<<tempLine.size()<<" templine="<<tempLine<<endl;
	 
    return(false);
  }else{
    return(false);
  }
}

//---------------------------------------------------------------------------
// * Copy
//---------------------------------------------------------------------------
//
void
Contig::CopyContig(Contig &NewC)
{
  Contig &C = (*this);

  cout<<"copy consensus...\t"<<flush;

  // Copy consensus. 
  NewC.resize(C.size());
  NewC.mName = C.mName;
  NewC.mDescription = C.mDescription;
  for(int c = 0;c< C.size();c++){
    NewC[c] = C[c];
  }
  
  cout<<"quality..\t"<<flush;

  // Copy quality
  NewC.mConsensusQuality.resize(C.mConsensusQuality.size());
  for(int c = 0;c< C.mConsensusQuality.size();c++){
    NewC.mConsensusQuality[c] = C.mConsensusQuality[c];
  }

  cout<<"data struc...\t"<<flush;

  // Copy data structures. 
  NewC.mPaddedBases = C.mPaddedBases;
  NewC.mNumEntries = C.mNumEntries;
  NewC.mNumBaseSeqs = C.mNumBaseSeqs;

  NewC.mTemplate2ReadMap = C.mTemplate2ReadMap;
  NewC.mRead2TemplateMap = C.mRead2TemplateMap;

  NewC.mReadTemplateDir = C.mReadTemplateDir;
  NewC.mReadIdxMap = C.mReadIdxMap;
  
  cout<<"reads...\t"<<flush;

  // Copy reads
  NewC.mReads.resize(C.mReads.size());
  for(int r = 0;r < C.mReads.size();r++){
    NewC.mReads[r] = C.mReads[r];
  }

  //  cout<<"readqual "<<flush;

  //cout<<C.mReadsQuality.size()<<"..\t"<<flush;

  //cout<<"("<<C.mReadsQuality[0].size()<<")"<<flush;

  // Copy reads quality. 
  NewC.mReadsQuality.resize(C.mReadsQuality.size());
  //  cout<<"("<<C.mReadsQuality[0].size()<<")"<<flush;
  for(int r = 0;r < C.mReadsQuality.size();r++){
    //cout<<C.mReadsQuality[r].size()<<" : "<<flush;
    NewC.mReadsQuality[r] = C.mReadsQuality[r];
  }


  //  cout<<"more info...\t"<<flush;

  NewC.mTInfo = C.mTInfo;
  NewC.mTemplateDir = C.mTemplateDir;
  NewC.mAlignStart = C.mAlignStart;
  NewC.mReadSense = C.mReadSense;
  NewC.mReadTrimLeft = C.mReadTrimLeft;
  NewC.mReadTrimRight = C.mReadTrimRight;
  NewC.mBaseSegments = C.mBaseSegments;
  NewC.mQALines = C.mQALines;
  NewC.mTrimLeft = C.mTrimLeft;
  NewC.mTrimRight = C.mTrimRight;
  NewC.mLeftmostBAC = C.mLeftmostBAC;
  NewC.mRightmostBAC = C.mRightmostBAC;

  NewC.mBACTrimLeft = C.mBACTrimLeft;
  NewC.mBACTrimRight = C.mBACTrimRight;
}


//----------------------------------------
// * CleanInitFromContigs
//----------------------------------------
// Recompute overall AceFile statistics and data structures from 
// the contigs.  Should be called when contigs have been added or 
// deleted to update AceFile info.  
// 
void 
AceFile::CleanInitFromContigs()
{
  BuildReadHash();

  ComputeAllContigStats();
}


//----------------------------------------
// * BuildReadHash
//----------------------------------------
// Creates a hash which maps read names to the indices of 
// the contigs which contain them. 
// 
void 
AceFile::BuildReadHash()
{
  mReadContigMap.clear();
	
  for(int c= 0;c< size();c++){
    Contig &C=(*this)[c];
    for(int r= 0;r< C.mReads.size();r++){
      mReadContigMap.Add(C.mReads[r].mName,c);
      C.mReadIdxMap[C.mReads[r].mName] = r;
    }
  }
}


//----------------------------------------
// * CalculateBACContigTrims
//----------------------------------------
// This ought to be good enough, but if it is ever in doubt we could read in a 
// list of BAC project names and compare these to that list. 
// 
//
bool 
AceFile::IsBACName(string Name){
  if (Name.size() < 10) return(false);

  // first five characters should be alpha
  for(int i = 0;i< 5; i++){
    if (!isalpha(Name[i])) return(false);
  }

  if (!isdigit(Name[5])) return(false);

  if (!isalpha(Name[6])) return(false);

  for(int i = 7;i< 11;i++){
    if (!isdigit(Name[i])) return(false);
  }

  // if all these conditions have been met, it's a BAC read. 
  return(true);
}

/*
//----------------------------------------
// * GetReadExtents
//----------------------------------------
// Given two contig positions, it finds all reads overlappign those two 
// contig positions and returns the leftmost and rightmost extent of any 
// of those reads. This is a utility for when we are considering merging
// contigs, and find the consensus overlap and want to 
//
// CStart,CEnd, leftmostRead,rightmostRead, are all given in depadded 
// coordinates.  
//
void
Contig::GetReadExtents(int CStart,int CEnd,int &leftmostRead,int &rightmostRead,
SequenceVector &Reads,SequenceQualityVector &Quals)
{

// CStart and CEnd are one based and in depadded coordinates.  Convert this to 
// padded zero based coordinates. 
//
// 1234556789011234
// 0123456789012345
// ACTG*GATCC*C*AAT
//              ^
// Depadded = 12
// Padded = 13  


leftmostRead = C.size();
rightmostRead = -1;

for(int r = 0;r < mReads.size();r++){
int rLeft = RIdx2CIdx(0,r);
int rRight = RIdx2CIdx(mReads[r].size()-1,r);

//
//	 CSt
// -------------           ----------
//	
//	           ---------     CEnd
//	
// Does this read overlap the CStart/CEnd region?
		
bool bOverlap = true;
if ((rLeft < CStart) && (rRight < CStart)) bOverlap = false;
if ((rLeft > CEnd) && (rRight > CEnd)) bOverlap = false;
	 
// If it overlaps, see if this is the rightmost or leftmost, 
// and save the read. 
if (bOverlap){
if (rLeft < leftmostRead) leftmostRead = rLeft;
if (rRight > rightmostRead) rightmostRead = rRight;	
Reads.push_back(mReads[r]);
Reads.push_back(mReadsQual[r]);
}
}
}
*/

//----------------------------------------
// * CalculateBACContigTrims
//----------------------------------------
// Finds the leftTrim and rightTrim values for each 
// contig where leftTrim and rightTrim are defined as the 
// most extreme left or right position in the contig that is
// overlapper by a BAC read.   That is, everything outside of 
// that region should be WGS only.   Relies on a function 
// IsBACName which is currently hard-coded for Rat projects. 
// 
void
AceFile::CalculateBACContigTrims(int phrapCutoff)
{
  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];
    C.mTrimLeft = -1;
    C.mTrimRight = -1;

    // Find left trim  (should replace with sliding window)
    for(int i = 0;i< C.size();i++){		
      if (mContigQualityVec[c][i] >= phrapCutoff) {
	C.mTrimLeft = i;
	break;
      }
    }

    // Find right trim (should replace with sliding window)
    for(int i = C.size()-1;i >=0;i--){
      if (mContigQualityVec[c][i]  >= phrapCutoff){
	C.mTrimRight = i;
	break;
      }
    }

    int leftmostBAC = C.size();  // impossible value to right
    int rightmostBAC = -1;                  // impossible value to left

    if ((C.mTrimLeft == -1) && (C.mTrimRight == -1)){
      C.mLeftmostBAC = -1;
      C.mRightmostBAC = -1;
    }else{
      for(int r = 0;r < C.mReads.size();r++){
	if (IsBACName(C.mReads[r].mName)){
	  // See if it's the leftmost BAC read
	  int leftConsensus = C.mReadTrimLeft[r] + C.mAlignStart[r]; 
	  if (leftConsensus >= 0){
	    if (leftConsensus < leftmostBAC) leftmostBAC = leftConsensus;
	  }
			 
	  // See if it's the rightmost BAC read.
	  int rightConsensus = C.mReadTrimRight[r] + C.mAlignStart[r]; 
	  if (rightConsensus <= (C.size() - 1)){
	    if (rightConsensus > rightmostBAC) rightmostBAC = rightConsensus;
	  }
	}
      }
		
      C.mLeftmostBAC = leftmostBAC;
      C.mRightmostBAC = rightmostBAC;
    }

    //cout<<C.mName<<"\tTrimLeft="<<C.mTrimLeft<<"\tTrimRight="<<C.mTrimRight;
    //cout<<"\tLeftmostBAC="<<C.mLeftmostBAC<<"\tRightmostBAC="<<C.mRightmostBAC<<endl;


    // If there is a valid BAC giving a leftmostBAC and rightmostBAC, 
    if ((C.mLeftmostBAC != -1) && (C.mRightmostBAC < C.size())){
      // determine if trim or BAC dominates
      C.mBACTrimLeft = max(C.mTrimLeft,C.mLeftmostBAC);
      C.mBACTrimRight = min(C.mTrimRight,C.mRightmostBAC);
    }else{
      // There is no valid BAC trim.  Maybe no BACs, but in any case once
      // the BACs are quality trimmed there is no BAC trim.  
      C.mBACTrimLeft = -1;
      C.mBACTrimRight = -1;
    }								  
  }
}


//----------------------------------------
// * CalculateTrim
//----------------------------------------
void
Contig::CalculateTrim(int cutoff){

  Contig &C = (*this);
  C.mTrimLeft = -1;
  C.mTrimRight = -1;

  // Find left trim
  for(int i = 0;i< C.size();i++){
    if (C.mConsensusQuality[i] >= cutoff) {
      C.mTrimLeft = i;
      break;
    }
  }
  
  // Find right trim
  for(int i = C.size()-1;i >=0;i--){
    if (C.mConsensusQuality[i]  >= cutoff){
      C.mTrimRight = i;
      break;
    }
  }	 
}

//----------------------------------------
// * CalculateContigTrims
//----------------------------------------
// Computes the left and right trim value for every contig, 
// where contigs are trimmed at the first base from the left
// or right where the consensus quality is >= cutoff
// 
void
AceFile::CalculateContigTrims(int cutoff)
{
  for(int c = 0;c< size();c++){	 
    Contig &C = (*this)[c];
    C.CalculateTrim(cutoff);
  }
}

//-----------------------------------------------------
// * CalculateReadTrims
//-----------------------------------------------------
// This is not a very good trim for assembly.  Some 
// sort of sliding window would be better.  For assessing 
// the extent of coverage compared to a poisson model, 
// though, this is probably about right.  
//
void
AceFile::CalculateReadTrims(int cutoff)
{

  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];

    // Sanity test.
    if (C.mReads.size() != C.mReadsQuality.size()){
      cerr<<C.mName<<"\tERROR: reads.size="<<C.mReads.size()<<" qualsize="<<C.mReadsQuality.size()<<endl;
      continue;
    }


    C.mReadTrimLeft.resize(C.mReads.size());
    C.mReadTrimRight.resize(C.mReads.size());
    for(int r = 0;r < C.mReads.size();r++){
      C.mReadTrimLeft[r] = -1;
      C.mReadTrimRight[r] = -1;

      if (C.mReads[r].size() != C.mReadsQuality[r].size()){
	cerr<<"ERROR: read "<<C.mReads[r].mName<<" reads.size="<<C.mReads[r].size()<<" qualsize="<<C.mReadsQuality[r].size()<<endl;
			 
	// If it's just empty, make up something and go on...
	if (C.mReadsQuality[r].size() == 0){
	  C.mReadsQuality[r].resize(C.mReads[r].size());
	  for(int i = 0;i< C.mReads[r].size();i++) C.mReadsQuality[r][i] = 20;
	}else{
	  // If the sizes mismatch, just skip this read. 
	  continue;
	}
      }

      // Find left trim
      for(int i = 0;i< C.mReads[r].size();i++){
	if (C.mReadsQuality[r][i] >= cutoff) {
	  C.mReadTrimLeft[r] = i;
	  break;
	}
      }
      // Find right trim
      for(int i = C.mReads[r].size()-1;i >=0;i--){
	if (C.mReadsQuality[r][i] >= cutoff){
	  C.mReadTrimRight[r] = i;
	  break;
	}
      }
		 
      //	cout<<C.mName<<" "<<C.mReads[r].mName<<"\t";
      //cout<<"ReadTrimLeft (qualonly)="<<C.mReadTrimLeft[r];
      //cout<<"ReadTrimRight (qualonly)="<<C.mReadTrimRight[r]<<endl;
		 
      int startX=-1,lastX=-1;
      if (C.mReads[r].size() > 50){
			 
	// See if there is an x at the start and/or end.   
	for(int i = 0;i< 50;i++) {
	  if (toupper(C.mReads[r][i]) == 'X'){
	    startX=i; 
	    break;
	  }
	}
			 
			 
	for(int i = C.mReads[r].size()-1;i > (C.mReads[r].size()-50);i--){
	  if (toupper(C.mReads[r][i]) == 'X'){
	    lastX=i; 
	    break;
	  }
	}
			 
	if (startX>=0){
	  for(int i = startX;i< C.mReads[r].size();i++){
	    if (toupper(C.mReads[r][i]) != 'X'){
	      // see if it's further in than the quality trim
	      if (i > C.mReadTrimLeft[r]){
		C.mReadTrimLeft[r] = i;
		break;
	      }
	    }
	  }
	}
			 
	if (lastX>=0){
	  for(int i = lastX;i>=0;i--){
	    if (toupper(C.mReads[r][i]) != 'X'){
	      // see if it's further in than the quality trim
	      if (i < C.mReadTrimRight[r]){
		C.mReadTrimRight[r] = i;
		break;
	      }
	    }
	  }
	}
      }
		 
      //cout<<C.mName<<" "<<C.mReads[r].mName<<"\t";
      //cout<<"ReadTrimLeft (final)="<<C.mReadTrimLeft[r];
      //cout<<"ReadTrimRight (final)="<<C.mReadTrimRight[r]<<endl;
		 
    } // for r
  } // for c
}

//----------------------------------------
// * ParseCO
//----------------------------------------
// Part of parser.  Parses the CO portion of the 
// AceFile, including the CO line itself and the 
// consensus sequence. 
// 
// Function assumes that line contains CO token. 
// and that fin points to the line after the CO line in
// the file. 
//
string
AceFile::ParseCO(istream &fin,string &line,int &ContigNum)
{

  // cout<<"ParseCO: \t"<<ContigNum<<" this.size="<<(*this).size()<<"\t";
  //cout<<line<<endl;

  if (ContigNum >= (*this).size()){
    cerr<<"ContigNum="<<ContigNum<<" (*this).size="<<(*this).size()<<endl;
    exit(1);
  }
  
  Contig &C = (*this)[ContigNum];
  C.clear();
  C.mReads.clear();
  C.mAlignStart.clear();
  C.mReadSense.clear();
  
  // ... summary contig info
  int offset = 2;
  C.mName = GetNextWSDelimitedString(line,offset,offset);
  C.mPaddedBases = GetNextInteger(line,offset,offset);
  C.mNumEntries = GetNextInteger(line,offset,offset);
  C.mNumBaseSeqs = GetNextInteger(line,offset,offset);
  
  // ... read consensus sequence...
  do{

    GetLine(fin,line);
    if (ContainsWSPadded(line,"BQ")) break;
    if (ContainsOnlyWS(line)) break;

    //	 StripWhiteSpace(line);	// KJD Too slow
    //  StripTrailingSpace(line);  // KJD faster but untested. 
    C.ConcatenateString(line);
    //C.ConcatenateNonWhiteSpace(line);
  }while(!fin.eof()); // Ends when BQ is found or EOF 
  
  // Small sanity test
  if (C.mPaddedBases != C.size()){	 
    cerr<<"WARNING: PaddedBases = "<<C.mPaddedBases;
    cerr<<" C.mConsensus.size="<<C.size()<<endl;
  }  
  return(line);
}

//----------------------------------------
// * ParseBQ
//----------------------------------------
// Function assumes that line contains BQ token and
// that fin points to the point right after the BQ 
// in the file.  
//
// Note that the Ace format does not include quality
// values for pads (*).  I want mConsensusQuality to 
// be the same size as the number of padded bases, so 
// I add quality values (e.g. 0) in for the pads here. 
// When an Ace file is written back out, these pad qualities
// are removed again. 
// 
string
AceFile::ParseBQ(istream &fin,string &line,int &ContigNum)
{
  Contig &C = (*this)[ContigNum];
  C.mConsensusQuality.clear();
  
  int seqIdx = 0;
  // ... read consensus sequence...
  do{
    GetLine(fin,line);
    if (ContainsWSPadded(line,"AF")) break;
    if (ContainsOnlyWS(line)) break;

    // Perhaps a series of GetNextInt calls would be faster...

    vector<int> temp;
    GetIntsFromLine(line,temp);
    for(int i = 0;i< temp.size();i++){
      if (C[seqIdx] == '*'){
	// There might be a run of *'s, if so, go through 
	// them adding zeros		 
	while(C[seqIdx] == '*'){
	  C.mConsensusQuality.push_back(0);
	  seqIdx++;
	};
	// Now we are at a non-* base. 		  
	C.mConsensusQuality.push_back(temp[i]);
	seqIdx++;
      }else{
	C.mConsensusQuality.push_back(temp[i]);
	seqIdx++;
      }
    }
  }while(!fin.eof()); // Ends when BQ is found or EOF 

  return(line);
}

//----------------------------------------
// * ParseAF
//----------------------------------------
// Parses the AF lines for a particular contig.  
// 
// Function assumes that line contains AF token.  
// 
string
AceFile::ParseAF(istream &fin,string &line,int &ContigNum)
{
  //cout<<"ParseAF line="<<line<<endl;
  Contig &C = (*this)[ContigNum];
  int ReadIdx = 0;
  do{
    C.mReads.resize(ReadIdx+1);
    C.mAlignStart.resize(ReadIdx+1);
    C.mReadSense.resize(ReadIdx+1);
    int offset = 2;
    string TI = GetNextWSDelimitedString(line,offset,offset);
    //cout<<"AF TI="<<TI<<endl;

    string Sense = GetNextWSDelimitedString(line,offset,offset);
    C.mReadSense[ReadIdx] = Sense[0];
    C.mAlignStart[ReadIdx] = GetNextInteger(line,offset,offset);
	 
    ReadIdx++;
    GetLine(fin,line);
    //cout<<"Parse AF line="<<line<<endl;
    if (!ContainsWSPadded(line,"AF")) break;
  }while(!fin.eof());		
  //cout<<"end ParseAF"<<endl;

  return(line);
}


//----------------------------------------
// * ParseBS
//----------------------------------------
// Saves the BS lines for a contig.  The BS lines are not
// really interpreted by this parser.   The write function
// for AceFile simply writes these lines back out verbatim, 
// except when the contig has been modified, in which case 
// it currently writes out one BS line only.  
//
// 
// Function assumes that line contains BS token.  
// 
string
AceFile::ParseBS(istream &fin,string &line,int &ContigNum)
{
  //cout<<"ParseAF line="<<line<<endl;
  Contig &C = (*this)[ContigNum];
  int SegIdx = 0;
  do{
    C.mBaseSegments.resize(SegIdx+1);
    int offset = 2;
    C.mBaseSegments[SegIdx].mStart = GetNextInteger(line,offset,offset);
    C.mBaseSegments[SegIdx].mEnd = GetNextInteger(line,offset,offset);
    C.mBaseSegments[SegIdx].mName = GetNextWSDelimitedString(line,offset,offset);
	 
    SegIdx++;
    GetLine(fin,line);
    if (!ContainsWSPadded(line,"BS")) break;
  }while(!fin.eof());		

  //cout<<C.mName<<" BaseSegments.size="<<C.mBaseSegments.size()<<endl;

  return(line);
}


//----------------------------------------
// * ParseRD
//----------------------------------------
// Function assumes that line contains RD token.  
// RD  ReadName padded_bases whole_read_infos read_tags
//
// Parses the RD lines and subsequence read sequence. 
// The read sequences are stored with each contig. 
//
string
AceFile::ParseRD(istream &fin,string &line,int &ContigNum)
{
  //cout<<"Parse RD Initial line="<<line<<endl;
  Contig &C = (*this)[ContigNum];
  int ReadIdx = 0;
  do{
    //cout<<"ParseRD do{ line="<<line<<endl;
    if (ContainsWSPadded(line,"RD")){
      C.mReads.resize(ReadIdx+1);
      int offset = 2;
      C.mReads[ReadIdx].mName = GetNextWSDelimitedString(line,offset,offset);
      //cout<<"ReadName="<<C.mReads[ReadIdx].mName<<endl;
      do{
	GetLine(fin,line);
	if (ContainsOnlyWS(line)) break;
	else{
	  //				StripWhiteSpace(line);
	  //    StripTrailingSpace(line);
	  C.mReads[ReadIdx].ConcatenateString(line);
	  //C.mReads[ReadIdx].ConcatenateNonWhiteSpace(line);
	}
      }while(!fin.eof());
      //cout<<"readsize="<<C.mReads[ReadIdx].size()<<endl;
    }else if (ContainsWSPadded(line,"DS")){
      int spaceIdx = line.find(" ");
      C.mReads[ReadIdx].mDescription = line.substr(spaceIdx+1,line.size() -spaceIdx-1);

      C.RecordTemplateInfo(ReadIdx);

      // Once find the DS line, are done with the read.   
      ReadIdx++;
      //cout<<"description="<<C.mReads[ReadIdx].mDescription<<endl;

    }else if (ContainsWSPadded(line,"QA")){
      C.mQALines.push_back(line);
    }else if (ContainsWSPadded(line,"CO")) break;

    GetLine(fin,line);

  }while(!fin.eof());

  //  cout<<"CO line? -->"<<line<<endl;

  //cout<<" End ParseRD"<<endl;

  return(line);
}


//----------------------------------------
// * Open
//----------------------------------------
// Function for Java SWIG.  
//
// Allows Java to open the file and get the number of 
// contigs.  
//
int 
AceFile::Open(string FName){
  mZin;
  mZin.open(FName.c_str());
  if (mZin.fail()){
    cerr<<"Unable to open "<<FName<<" as an ace file."<<endl;
    exit(1);
  }

  string line;
  int offset = 0;
  GetLine(mZin,line);	 
  
  // We're at the beginning of the file.   
  if (line.find("AS") >=0){
    mNumContigs = GetNextInteger(line,offset,offset);
    mNumReads = GetNextInteger(line,offset,offset);
    resize(mNumContigs);
    mContigNum = 0;
    return(mNumReads);
  }else{
    cerr<<"ERROR: File does not start with AS. "<<endl;
    exit(1);
  }
}

//----------------------------------------
// * Close
//----------------------------------------
// Function for Java SWIG.  
//
// Allows Java to close the opened file stream. 
//
void
AceFile::Close(){
  mZin.close();
}

//----------------------------------------
// * ReadNextContig
//----------------------------------------
// Function for Java SWIG.  
//
// Reads the Ace file one contig at a time, to allow the
// calling Java program a chance to update the progress bar. 
//  
// KJD Normally I would not make mLine a member variable
// but would pass it in as an reference argument, but I want
// to avoid exposing that sort of thing to Java, so I hide
// it by making it a member variable. 
// 


//This should be ReadContigBatch()
//It should read contigs in batches of 50Kbp. 
//However, how can we tell Java how big it is in advance?  
//Maybe still do that based on # reads, but reads will jump by 
//50Kbp batch sizes. 

int 
AceFile::ReadContigBatch()
{
  int totalReads = 0;
  int numReads;
  do{
    if (FindTokenPlusWS(mZin,"CO",mLine)){
      mLine = ParseCO(mZin,mLine,mContigNum);
		  
      // BQ is optional, so if it is there, parse it, if not, go on. 
      if (FindTokenPlusWS(mZin,"BQ",mLine)){
	mLine = ParseBQ(mZin,mLine,mContigNum);
		  
	if (FindTokenPlusWS(mZin,"AF",mLine)){	 			
	  mLine=ParseAF(mZin,mLine,mContigNum);

	  // KJD This is too brittle.   Don't want to insist that BS is there. 
	  if (FindTokenPlusWS(mZin,"BS",mLine)){
	    mLine = ParseBS(mZin,mLine,mContigNum);
				
	    if (FindTokenPlusWS(mZin,"RD",mLine)){
	      mLine = ParseRD(mZin,mLine,mContigNum);
	    }
	  }
	}
      }		  
      mContigNum++;		
    }  
  
    numReads = (*this)[mContigNum-1].mReads.size();
    totalReads +=numReads;
  }while(!mZin.eof() && (totalReads < 250));

  if (mZin.eof()) return(-1);
  else return(totalReads);
}


//----------------------------------------
// * Read
//----------------------------------------
// Reads the Ace file proper.  
// 
void 
AceFile::Read(istream &fin){
  string line;
  int offset = 0;

  mbIsConsistent = true;

  GetLine(fin,line);	 
  if (line.find("AS") >=0){
    mNumContigs = GetNextInteger(line,offset,offset);
    //cout<<"NumContigs="<<mNumContigs<<endl;
    mNumReads = GetNextInteger(line,offset,offset);
    //cout<<"NumReads="<<mNumReads<<endl;
    resize(mNumContigs);
    int ContigNum = 0;
    string line;

    do{
      if (FindTokenPlusWS(fin,"CO",line)){
	line = ParseCO(fin,line,ContigNum);
		  
	// BQ is optional, so if it is there, parse it, if not, go on. 
	if (FindTokenPlusWS(fin,"BQ",line)){
	  line = ParseBQ(fin,line,ContigNum);

	  if (FindTokenPlusWS(fin,"AF",line)){	 			
	    line=ParseAF(fin,line,ContigNum);

	    if (FindTokenPlusWS(fin,"BS",line)){
	      line = ParseBS(fin,line,ContigNum);
			 
	      if (FindTokenPlusWS(fin,"RD",line)){
		line = ParseRD(fin,line,ContigNum);
	      }
	    }
	  }
	}		  
	//		  cout<<".";
	ContigNum++;		
      }			
    }while(!fin.eof());
  }else{
    cerr<<"File doesn't start with AS"<<endl;
    return;
  } 

  // Consistency check....
  int numContigs = 0;
  for(int i = 0;i< size();i++){
    Contig &C = (*this)[i];
    if (C.size() > 0){
      numContigs++;
    }

    if (C.mReads.size() != C.mNumEntries){
      cerr<<"INCONSISTENT CONTIG READS: "<<C.mName<<" Reads seen="<<C.mReads.size();
      cerr<<"  reads reported="<<C.mNumEntries<<endl;
      mbIsConsistent = false;
    }

  }
	
  if (numContigs != mNumContigs){
    cerr<<"INCONSISTENT ACE FILE: Contigs seen = "<<numContigs<<" contigs reported="<<mNumContigs<<endl;
    mbIsConsistent = false;
  }
}

void 
AceFile::Read(string FName,SequenceQualityVector &ReadQuals)
{
  izipstream zin;
  zin.open(FName.c_str());
  if (zin.fail()){
    cerr<<"Unable to open "<<FName<<" as an ace file."<<endl;
    return;
  }
  Read(zin);
  zin.close();
	
  BuildReadHash();
	
  // Time on big test ace (8.2MB compressed). without these following bits is 
  // Without = 11.7 - 12.0 sec. 
  // With =  13.3 - 13.5 sec.   (unzipped = 12.0 - 16.1 sec.).  
  CopyBQToQuality();

  AddReadsQuality(ReadQuals);
  
  cout<<"CalculateContigTrims..."<<flush;
  CalculateContigTrims(mContigTrim);
  cout<<"done."<<endl;
  cout<<"CalculateReadTrims..."<<flush;
  CalculateReadTrims(mPhredTrim);
  cout<<"done."<<endl;
  	
  ComputeAllContigStats();	
}

void 
AceFile::Read(string FName)
{
  izipstream zin;
  zin.open(FName.c_str());
  if (zin.fail()){
    cerr<<"Unable to open "<<FName<<" as an ace file."<<endl;
    return;
  }
  Read(zin);
  zin.close();

  BuildReadHash();

  // Time on big test ace (8.2MB compressed). without these following bits is 
  // Without = 11.7 - 12.0 sec. 
  // With =  13.3 - 13.5 sec.   (unzipped = 12.0 - 16.1 sec.).  
  CopyBQToQuality();

  // Goes through all of the reads and makes sure that they have some 
  // defined quality values, zeros are added where currently undefined. 
  // The purpose is to ensure that other code that expects the quality 
  // values to be there won't crash trying to access nonexistent values
  // (and it's easier to do this once than to test each instance to see if 
  // the quality exists). 
  AddFakeQualityValues();

  CalculateContigTrims(mContigTrim);
  CalculateReadTrims(mPhredTrim);

  ComputeAllContigStats();
	
}


//---------------------------------------------------
// * FindNewestAce
//---------------------------------------------------
// Looks for ace files in format AceRoot.k and 
// returns the full name of the one with the highest
// value for k. 
// 
bool 
AceFile::FindNewestAce(string &AceRoot)
{
  string cmd;

  cmd="rm -f /var/tmp/kjd-ace-stat-temp";
  system(cmd.c_str());

  cmd="ls ";
  cmd+=AceRoot;
  cmd+="* > /var/tmp/kjd-ace-stat-temp";
  system(cmd.c_str());

  vector<string> AceNames;
  string result ="/var/tmp/kjd-ace-stat-temp";  
  ReadFOF(result,AceNames);

  if ((AceNames.size() == 0) || (Contains(AceNames[0],"No such file"))){ 

  }

  sort(AceNames.begin(),AceNames.end());
  AceRoot = AceNames[AceNames.size()-1];

  if (Contains(AceRoot,".ace")) return(true);
  else return(false);
}


//--------------------------------------------------
// * GetNextQuality
//--------------------------------------------------
// Reads a single read quality from a multiple fasta quality file. 
//
bool
AceFile::GetNextQuality(istream &fin,SequenceQuality &Qual)
{
  Qual.mName ="";
  Qual.mDescription="";
  Qual.clear();

  // Must be first line  
  if (readbuf[0] == '^'){
    fin.getline(readbuf,READ_BUFFER_SIZE);
  }

  string temp = readbuf;
  
  // Extract name and description from line.  
  string nextDesc = readbuf;
  int offset = 1; // To skip >
  // Everything up to first white space is name, after that is description. 
  Qual.mName = GetNextWSDelimitedString(nextDesc,offset,offset);  

  // It's possible there is no description, so check:
  if ((offset >= 0) && (offset < nextDesc.size())){
    Qual.mDescription = nextDesc.substr(offset,nextDesc.size()-offset);
  }else{
    Qual.mDescription = "";
  }
  
  // Get the sequence until another > is encountered. 
  do{
    fin.getline(readbuf,READ_BUFFER_SIZE);	
    if (readbuf[0] == '>') break;
    GetIntsFromLine(readbuf,Qual);
  }while(!fin.eof());
  
  if (Qual.size()>0) return(true);
  return (false);
}


//---------------------------------------------------
// * AddReadsQuality
//---------------------------------------------------
// Reads the quality values for the reads out of a multiple
// fasta type quality file.   The function handles finding
// the contig and associated read for the quality, inserting
// pads, and reversing quality for reversed reads.  The 
// Ace file must be read first, before calling this function. 
//
void
AceFile::AddReadsQuality(SequenceQualityVector &Quals)
{	
  // Make all of the qualities match the size of the sequences
  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];
    C.mReadsQuality.resize(C.mReads.size());
  }

  cout<<"AceFile::Add "<<Quals.size()<<" read qualities."<<endl;
  int totalMissing = 0;
	
  for(int i = 0;i < Quals.size();i++){
		
    //		if (Quals[i].mName == "gnl|ti|1000001066"){
    //	cout<<"Well, gnl|ti|1000001066 is in the Quals input."<<endl;
    //}
		
    SequenceQuality &temp = Quals[i];

    // There might be several contigs with this read, add it to all 
    // that apply...
    vector<int> contigNums;
    mReadContigMap.GetValues(temp.mName,contigNums);

    for(int c = 0;c< contigNums.size();c++){
      int contigNum = contigNums[c];

      //if (Quals[i].mName == "gnl|ti|1000001066") cout<<"gnl|ti|1000001066 contigNum="<<contigNum<<endl;

      // Let's double check that these values actually exist
      if (contigNum < size()){
	Contig &C = (*this)[contigNum];				
	int readNum = C.mReadIdxMap[temp.mName];

	//if (Quals[i].mName == "gnl|ti|1000001066"){
	//	cout<<"gnl|ti|1000001066 readNum="<<readNum<<endl;
	//	}

	if (readNum < C.mReads.size()){			
	  // KJD If the size doesn't equal zero, we must have already 
	  // KJD read in this quality.   Something more robust is in order. 
	  if (C.mReadsQuality[readNum].size() == 0){
	    C.mReadsQuality[readNum].mName = temp.mName;
	    C.mReadsQuality[readNum].mDescription = temp.mDescription;
						
	    // Reverse the quality if the read is complemented
	    if (C.mReadSense[readNum] == 'C') temp.Reverse();
						
	    int idx = 0;
	    for(int s = 0;s< C.mReads[readNum].size();s++){
	      if (C.mReads[readNum][s] =='*') {
		//if (Quals[i].mName == "gnl|ti|1000001066") cout<<"0 ";
		C.mReadsQuality[readNum].push_back(0);
	      }else {
		if (idx< temp.size()){
		  //if (Quals[i].mName == "gnl|ti|1000001066") cout<<temp[idx]<<" ";
		  C.mReadsQuality[readNum].push_back(temp[idx++]);
		}else{
		  cerr<<C.mReads[readNum].mName<<" Hey tempsize="<<temp.size();
		  cerr<<" idx="<<idx<<"\t readsize=";
		  cerr<<C.mReads[readNum].size()<<endl;
		}
	      }
	    }
	  }else{
	    //						if (Quals[i].mName == "gnl|ti|1000001066") cout<<"Size NOT zero.  Quals[i].mName.size="<<Quals[i].mName.size()<<endl;
	  }
	}
	//				if (Quals[i].mName == "gnl|ti|1000001066") cout<<"C.mReads[readNum].size="<<C.mReads[readNum].size()<<endl;
      }
    } // for each contig containing the read.
  } // for each input quality...	
}




//---------------------------------------------------
// * GetReadsQuality
//---------------------------------------------------
// Reads the quality values for the reads out of a multiple
// fasta type quality file.   The function handles finding
// the contig and associated read for the quality, inserting
// pads, and reversing quality for reversed reads.  The 
// Ace file must be read first, before calling this function. 
//
void
AceFile::GetReadsQuality(istream &fin)
{
  // Make all of the qualities match the size of the sequences
  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];
    C.mReadsQuality.resize(C.mReads.size());
  }
  
  readbuf[0]='^';  // Indicates first read to GetNextQuality
  SequenceQuality temp;
  do{
    if (GetNextQuality(fin,temp)){
      // Now find where this one goes, if anywhere, and copy the 
      // contig values in with pads.  It might go several places...
      vector<int> contigNums;
      mReadContigMap.GetValues(temp.mName,contigNums);
      for(int i = 0;i < contigNums.size();i++){
	int contigNum = contigNums[i];

	// Let's double check that these values actually exist
	if (contigNum < size()){
	  Contig &C = (*this)[contigNum];
					
	  int readNum = C.mReadIdxMap[temp.mName];
					
	  if (readNum < C.mReads.size()){			
						
	    // KJD If the size doesn't equal zero, we must have already 
	    // KJD read in this quality.   Something more robust is in order. 
	    if (C.mReadsQuality[readNum].size() == 0){
	      C.mReadsQuality[readNum].mName = temp.mName;
	      C.mReadsQuality[readNum].mDescription = temp.mDescription;
							
	      // Reverse the quality if the read is complemented
	      if (C.mReadSense[readNum] == 'C') temp.Reverse();
							
	      int idx = 0;
	      for(int s = 0;s< C.mReads[readNum].size();s++){
		if (C.mReads[readNum][s] =='*') C.mReadsQuality[readNum].push_back(0);
		else {
		  if (idx< temp.size()){
		    C.mReadsQuality[readNum].push_back(temp[idx++]);
		  }else{
		    cerr<<C.mReads[readNum].mName<<" Hey tempsize="<<temp.size();
		    cerr<<" idx="<<idx<<"\t readsize=";
		    cerr<<C.mReads[readNum].size()<<endl;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }while(!fin.eof()); 
}



//----------------------------------------
// * gzExists
//----------------------------------------
bool
AceFile::gzExists(string FileName)
{
  string gzname = FileName+".gz";
  ifstream testin;
  testin.open(gzname.c_str());
  if (testin.fail()) return(false);
  else return(true);
}


//----------------------------------------
// * ReadAll
//----------------------------------------
// Reads Acefile and associated contig and read quality files. 
// Derives names from a root name. 
// 
void
AceFile::ReadAll(string &AceRoot) 
{
  FASTAFile FA;

  // Read the stuff that is in the ace file.  Basically 
  // this is everything except the read quality.  

  string AceName = AceRoot+".ace";
  string ContigQualName = AceRoot+".contigs.qual";
  string ReadQualName = AceRoot+".qual";  
  
  //cerr<<"Reading ";
  if (gzExists(AceName)){
    izipstream zin;
    AceName+=".gz";
    //cerr<<AceName<<"...";
    zin.open(AceName.c_str());
    Read(zin);
    zin.close();
  }else{
    ifstream fin;
    fin.open(AceName.c_str());
    //cerr<<AceName<<"...";
    Read(fin);
    fin.close();
  }
  //cerr<<"done."<<endl;

  // Now I need a hash that maps read names to the contig and read number
  // within that contig.  
  BuildReadHash();

  //cerr<<"Reading ";
  if (gzExists(ReadQualName)){
    izipstream zin;
    ReadQualName+=".gz";
    //cerr<<ReadQualName<<"...";
    zin.open(ReadQualName.c_str());
    GetReadsQuality(zin);
    zin.close();
  }else{
    ifstream fin;
    fin.open(ReadQualName.c_str());
    if (!fin.fail()){
      //cerr<<ReadQualName<<"...";
      GetReadsQuality(fin);
      fin.close();
    }
  }
  //  cerr<<"done."<<endl;

  CopyBQToQuality();

  // Goes through all of the reads and makes sure that they have some 
  // defined quality values, zeros are added where currently undefined. 
  // The purpose is to ensure that other code that expects the quality 
  // values to be there won't crash trying to access nonexistent values
  // (and it's easier to do this once than to test each instance to see if 
  // the quality exists). 
  AddFakeQualityValues();

  CalculateReadTrims(mPhredTrim);
  CalculateContigTrims(mContigTrim);  
}


//----------------------------------------
// * ReadAll
//----------------------------------------
// Reads Acefile and associated contig and read quality files. 
// Three file names given explicitly. 
// 
void
AceFile::ReadAll(string &AceName,string &ReadQualName) 
{
  FASTAFile FA;
	
  cerr<<endl;
  cerr<<"Reading "<<AceName<<"..."<<flush;
  if (Contains(AceName,".gz")){
    izipstream zin;
    zin.open(AceName.c_str());
    Read(zin);
    zin.close();
  }else{
    ifstream fin;
    fin.open(AceName.c_str());
    Read(fin);
    fin.close();
  }
  cerr<<"done."<<endl;

  // Now I need a hash that maps read names to the contig and read number
  // within that contig.  
  BuildReadHash();

  cerr<<"Reading "<<ReadQualName<<"..."<<flush;
  if (Contains(ReadQualName,".gz")){
    izipstream zin;
    zin.open(ReadQualName.c_str());
    GetReadsQuality(zin);
    zin.close();
  }else{
    ifstream fin;
    fin.open(ReadQualName.c_str());
    GetReadsQuality(fin);
    fin.close();
  }
  cerr<<"done."<<endl;

  CopyBQToQuality();

  // Goes through all of the reads and makes sure that they have some 
  // defined quality values, zeros are added where currently undefined. 
  // The purpose is to ensure that other code that expects the quality 
  // values to be there won't crash trying to access nonexistent values
  // (and it's easier to do this once than to test each instance to see if 
  // the quality exists). 
  AddFakeQualityValues();

  CalculateReadTrims(mPhredTrim);
  CalculateContigTrims(mContigTrim);
}


//----------------------------------------
// * CopyBQToQuality
//----------------------------------------
// This is a temporary hack as a bridge between the old
// way of getting the quality, from contigs.qual file, and 
// the new way, directly from the ACE file. 
//
// KJD DEPRICATE
//
void
AceFile::CopyBQToQuality()
{
  // Temporary hack. 

  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];

    // If we've read one already
    if (mContigQualityVec.size() > c){
      mContigQualityVec[c].clear();
    }else{
      mContigQualityVec.resize(mContigQualityVec.size()+1);
    }
    mContigQualityVec[c].resize(C.mConsensusQuality.size());
    for(int i = 0;i< C.mConsensusQuality.size();i++){
      mContigQualityVec[c][i] = C.mConsensusQuality[i];
    }

    //cout<<C.mName<<" mContigQualityVec="<<mContigQualityVec[c].size()<<" ConsensusQual=";
    //cout<<C.mConsensusQuality.size()<<endl;
  }
}


//----------------------------------------
// * AddFakeQuality
//----------------------------------------
// The data is never perfect, is it?  Basically this function ensures
// that the size of the quality data structure matchaes that of the 
// read, filling in zero quality if needed.   
// 
void 
AceFile::AddFakeQualityValues()
{
  for(int c = 0;c< size();c++){
    Contig &C = (*this)[c];

    // First, fixup the contig quality if needed.  Sadly, even this doesn't come
    // out right sometimes. 

    if (C.size() != mContigQualityVec[c].size()){
      // Don't want to bug people if there was no read quality at all. 
      if (mContigQualityVec[c].size() != 0){
	cerr<<"Forcing quality size match:";
	cerr<<C.mName<<" size="<<C.size()<<" qual size="<<mContigQualityVec[c].size();
	cerr<<" BQ.size="<<C.mConsensusQuality.size()<<endl;
      }
      mContigQualityVec[c].resize(C.size());
      for(int j = 0;j< mContigQualityVec[c].size();j++){
	mContigQualityVec[c][j] = 20;  // Some arbitrary quality designed not to scuttle the whole contig. 
      }
    }

    if (C.mReadsQuality.size() != C.mReads.size()){
      if (C.mReadsQuality.size() != 0){
	cerr<<"Forcing size match: ";
	cerr<<C.mName<<" read size="<<C.mReads.size();
	cerr<<C.mName<<" quality size="<<C.mReadsQuality.size()<<endl;
      }
      C.mReadsQuality.resize(C.mReads.size());
    }

    for(int r = 0;r < C.mReads.size();r++){
      if (C.mReads[r].size() != C.mReadsQuality[r].size()){
	// No need to bug if there is no quality to begin with. 
	if (C.mReadsQuality[r].size() != 0){
	  cerr<<"\t Forcing size to match: ";
	  cerr<<C.mReads[r].mName<<" readsize="<<C.mReads[r].size()<<" qual size="<<C.mReadsQuality[r].size()<<endl;
	}
	C.mReadsQuality[r].resize(C.mReads[r].size());
	for(int i = 0;i< C.mReads[r].size();i++){
	  C.mReadsQuality[r][i] = 0;
	}
      }
    }
  }
}





//---------------------------------------------------
// * ComputeAllcontigStats
//---------------------------------------------------
// Computes basic statistics for the trimmed portion of 
// a contig. 
// 
void 
AceFile::ComputeAllContigStats()
{
  mTotalContigLength = 0;
  mContigsMeetingCriteria = 0;
  mSizeMeetingCriteria = 0;

  mTrimmedContigLength = 0;
  mWGSReadCount = 0;
  mTotalReadCount = 0;

  mMaxContigSize = 0;
  mMaxTrimmedContigSize = 0;

  // Go through each contig, computing some statistics. 
  for(int i = 0;i< size();i++){
    Contig &C = (*this)[i];
    int contigSize = C.size();
    if (contigSize > mMaxContigSize) mMaxContigSize = contigSize;
    int trimmedSize = (C.mTrimRight-C.mTrimLeft+1);
    if (trimmedSize > mMaxTrimmedContigSize) mMaxTrimmedContigSize = trimmedSize;
    mTotalContigLength +=contigSize;
    mTrimmedContigLength += trimmedSize;

    if (contigSize >=mMinContigLength){
      mContigsMeetingCriteria++;
      mSizeMeetingCriteria+=contigSize;

      for(int r = 0;r< C.mReads.size();r++){
	int idx = C.mReads[r].mName.find("gnl|ti|");
	if (idx >=0){
	  mWGSReadCount++;
	}
      }
      mTotalReadCount+=C.mReads.size();		
    }
	 
    if (trimmedSize >=mMinContigLength){
      mTrimmedContigsMeetingCriteria++;
      mTrimmedSizeMeetingCriteria+=trimmedSize;
    }
  }

  mN50 = N50();
  mN50Trimmed = N50Trimmed(); 

  CalcAverageCoverage();
  CalcAverageCoverageTrimmed();
}

//---------------------------------------------------
// * N50Trimmed
//---------------------------------------------------
// Computes the N50 Size for a trimmed contig. 
// 
int 
AceFile::N50Trimmed()
{
  vector<int> sizeList;
  
  GetSortedSizesTrimmed(sizeList);
 
  int totalBases = mTrimmedContigLength;
  int halfBases = (int)((float)totalBases /(float) 2.0);
  int sizeSum = 0;
  for(int i = 0;i< sizeList.size();i++){	
    sizeSum += sizeList[i];
    if (sizeSum >= halfBases) return(sizeList[i]);
  }
}

//---------------------------------------------------
// * GetSortedSizesTrimmed
//---------------------------------------------------
//  Helper function for N50Trimmed
void 
AceFile::GetSortedSizesTrimmed(vector<int> &sizeList)
{
  for(int i = 0;i< size();i++){
    Contig &C = (*this)[i];
    int trimmedSize = (C.mTrimRight-C.mTrimLeft+1);
    if (trimmedSize >= mMinContigLength){
      sizeList.push_back(trimmedSize);
    }
  }
  sort(sizeList.begin(),sizeList.end(),Ace_lt_int);
}

//---------------------------------------------------
// * N50
//---------------------------------------------------
// Computes the N50 size for untrimmed contigs. 
// 
int 
AceFile::N50()
{
  vector<int> sizeList;
  
  GetSortedSizes(sizeList);
  
  int totalBases = mTotalContigLength;
  int halfBases = (int)((float)totalBases /(float) 2.0);
  int sizeSum = 0;
  for(int i = 0;i< sizeList.size();i++){	
    sizeSum += sizeList[i];
    if (sizeSum >= halfBases) return(sizeList[i]);
  }
}

//---------------------------------------------------
// * GetSortedsizes
//---------------------------------------------------
//
void 
AceFile::GetSortedSizes(vector<int> &sizeList)
{
  for(int i = 0;i< size();i++){
    if ((*this)[i].size() >= mMinContigLength){
      sizeList.push_back((*this)[i].size());
    }
  }
  sort(sizeList.begin(),sizeList.end(),Ace_lt_int);
}

//----------------------------------------
// * CalcAverageCoverage()
//----------------------------------------
//
void
AceFile::CalcAverageCoverage()
{
  int totalCoverage = 0;
  int totalSize = 0;
  for(int i = 0;i< size();i++){
    if ((*this)[i].size() >= mMinContigLength){
      totalSize +=(*this)[i].size();
      for(int c = 0;c< (*this)[i].size();c++){
	int coverage = (*this)[i].CalcCoverageTrimmed(c);
	totalCoverage+=coverage;
      }
    }
  }
  mAverageCoverage = (float)totalCoverage/(float)totalSize;
}

//----------------------------------------
// * CalcAverageCoverageTrimmed()
//----------------------------------------
//
void
AceFile::CalcAverageCoverageTrimmed()
{
  int totalCoverage = 0;
  int totalSize = 0;
  for(int i = 0;i< size();i++){
    int trimmedSize = (*this)[i].mTrimRight - (*this)[i].mTrimLeft+1;
    if (trimmedSize >= mMinContigLength){
      totalSize +=trimmedSize;
      for(int c = (*this)[i].mTrimLeft;c<= (*this)[i].mTrimRight;c++){		  
	int coverage = (*this)[i].CalcCoverageTrimmed(c);
	totalCoverage+=coverage;
      }
    }
  }
  mTrimmedAverageCoverage = (float)totalCoverage/(float)totalSize;
}


//----------------------------------------
// * PrefixContigs
//----------------------------------------
// Adds some prefix to the contig names. 
// Used when combining ace files, typically 
// to produce a name like:
//
// proj_gaab_Contig1
// 
void
AceFile::PrefixContigs(string prefix)
{
  for(int i = 0;i< size();i++){
    (*this)[i].mName = prefix+(*this)[i].mName;
  }
}

void
AceFile::RemoveContig(int c)
{
  int removedReads = (*this)[c].mReads.size();

  (*this).erase((*this).begin()+c);
  mContigQualityVec.erase(mContigQualityVec.begin()+c);
  
  mNumReads -=removedReads;

  mNumContigs--;
}

//----------------------------------------
// * GetContigsSharingTemplates
//----------------------------------------
void 
AceFile::GetContigsSharingTemplates(int cIdx,vector<int> &ConnectedContigs)
{

  // Put them into a set to get a list with no duplicates. 
  IntSet ConnectedSet;

  StringSet sourceTID;
  Contig &C = (*this)[cIdx];

  int bugcsize = C.size();


  // Make a set of the source template IDs
  C.GetTemplateInfo(C.mTInfo);
  for(int i = 0;i< C.mTInfo.size();i++){
    sourceTID.Add(C.mTInfo[i].mTID);
  }

  // Make sure we have template info for everything. 
  for(int c = 0;c< size();c++){
    if (c != cIdx){
      Contig &TestC = (*this)[c];		
      TestC.GetTemplateInfo(TestC.mTInfo);
    }
  }


  for(int c = 0;c< size();c++){
    if (c != cIdx){
      Contig &TestC = (*this)[c];
      for(int i = 0;i< TestC.mTInfo.size();i++){
	if (sourceTID.Contains(TestC.mTInfo[i].mTID)){
	  ConnectedSet.Add(c);
	}
      }
    }
  }

  IntSet::iterator it;
  for(it = ConnectedSet.begin(); it != ConnectedSet.end();it++){
    ConnectedContigs.push_back((*it));
  }
}

//----------------------------------------
// * FindTemplatesBetween
//----------------------------------------
// 
void 
AceFile::FindTemplatesBetween(Contig &C1, Contig &C2, vector<TemplateInfo> &TInfo)
{
  // Get the template info from the read descriptions
  C1.GetTemplateInfo(C1.mTInfo);
  C2.GetTemplateInfo(C2.mTInfo);
	
  // Create sets of template ID for the two contigs
  StringSet C1Templates,C2Templates;
  for(int i = 0;i< C1.mTInfo.size();i++){
    C1Templates.Add(C1.mTInfo[i].mTID);
  }
  for(int i = 0;i< C2.mTInfo.size();i++){
    C2Templates.Add(C2.mTInfo[i].mTID);
  }
	
  // Find out which templates occurr in both contigs
  StringSet BothTemplates;
  StringSet::iterator it;
  for(it = C1Templates.begin();it != C1Templates.end();it++){
    string tid = (*it);
    if (C2Templates.Contains(tid)){
      BothTemplates.Add(tid);
    }
  }
  
  // Save template info for templates found in both. 
  // Arbitrarialy return the one from c1. 
  for(it = BothTemplates.begin();it != BothTemplates.end();it++){
    string tid = (*it);
    for(int t = 0;t< C1.mTInfo.size();t++){
      if (Contains(C1.mTInfo[t].mTID,tid)){
	TInfo.push_back(C1.mTInfo[t]);
      }
    }
  }	
}


//----------------------------------------
// * FindTemplatesBetween
//----------------------------------------
// Finds all of the templates that connect contig c1 with contig c2 (if any). 
// The TemplateInfo entry returned for each pair is that of the c1 read. 
// 
void 
AceFile::FindTemplatesBetween(int c1,int c2,vector<TemplateInfo> &TInfo)
{
  if ((c1 < 0) || (c1 >= size())) return;
  if ((c2 < 0) || (c2 >= size())) return;

  Contig &C1 = (*this)[c1];
  Contig &C2 = (*this)[c2];

  FindTemplatesBetween(C1,C2,TInfo);
}


//=======================================================================================
// --------------------------------   Contig  -------------------------------------------
//=======================================================================================


//----------------------------------------
// * CountReadsMatching
//----------------------------------------
// Counts how many read names in this contig match the key. 
// Eventually this will be changed to a regular expression.
// Use, for example, to see how many wgs reads there are 
// in a contig. 
// 
int
Contig::CountReadsMatching(string key)
{
  int count = 0;
  for(int r = 0;r< mReads.size();r++){
    if (Contains(mReads[r].mName,key)) count++;
  }
  return(count);
}


//----------------------------------------
// * AllReadsMatch
//----------------------------------------
// Determines if all of the reads 
// 
bool
Contig::AllReadsMatch(string key)
{
  int numReadsMatched = CountReadsMatching(key);
  if (numReadsMatched == mReads.size()) return(true);
  else return(false);
}



//----------------------------------------
// * NoReadsMatch
//----------------------------------------
// Determines if none of the reads match. 
// 
bool
Contig::NoReadsMatch(string key)
{
  int numReadsMatched = CountReadsMatching(key);
  if (numReadsMatched == 0) return(true);
  else return(false);
}


//----------------------------------------
// * ComputeSyntheticBS
//----------------------------------------
// When I compute a consensus I do it column by 
// column. So, strictly speaking, I do not assemble
// a consensus from segments of reads like Phrap does. 
// Still, all of the Phrap related tools may expect to 
// see BS lines that can be interpreted this way.  
// This function generates BS lines that are consistent
// with the reads and the consensus. 
// 
// 
void 
Contig::ComputeSyntheticBS()
{
  int nextBS = 0;

  int errcnt = 0;
  vector<int> examples;

  mBaseSegments.clear(); // Make sure there aren't any to start with. 
  mNumBaseSeqs = 0;

  // For each consensus position...
  for(int c = 0; c < size();c++){

    bool bFound = false;

    // Find a read that maps to that consensus position...
    for(int r = 0; r < mReads.size(); r++){
      int rIdx =  CIdx2RIdx(c,r);

      if ((rIdx >= 0) && (rIdx < mReads[r].size())){
	bFound = true;

	// If the read matches the consensus, or the consensus is a *,  start a run for a BS entry...
	if (toupper(mReads[r][rIdx]) == toupper((*this)[c]) ||
	    ((*this)[c] == '*') || 
	    (toupper((*this)[c]) == 'N') || 
	    (toupper((*this)[c]) == 'X')){

	  mBaseSegments.resize(nextBS+1);
	  mBaseSegments[nextBS].mStart = c+1;
	  mBaseSegments[nextBS].mName = mReads[r].mName;
					
	  // Continue the run for as long as the read matches...
	  while(toupper(mReads[r][rIdx]) == toupper((*this)[c]) || 
		((*this)[c] == '*') || 
		(toupper((*this)[c]) == 'N') || 
		(toupper((*this)[c]) == 'X')){
	    mBaseSegments[nextBS].mEnd = c+1;
	    rIdx++;
	    c++;
	    // Stop run if we are at end of consensus or end of the read...
	    if (c >= size()) break;
	    if (rIdx >= mReads[r].size()) break;
	  }
	  nextBS++;
	}
      }
    }
		
    // None of the reads seem to overlap this consensus position. 
    if (!bFound){
      errcnt++;
		  
      if (errcnt < 10){
	examples.push_back(c);
      }
    }
  }


  if (errcnt > 0){
    cerr<<"\n";
    cerr<<mName<<"\tSyntheticBS WARNING: Unable to find read for ";
    cerr<<errcnt<<" consensus positions."<<endl; 
    cerr<<"For example: "<<endl;
    for(int i = 0;i < examples.size();i++) cout<<"\tconsensus "<<examples[i]<<endl;

    cerr<<"mReads.size="<<mReads.size()<<endl;
    cerr<<"mAlignStart.size="<<mAlignStart.size()<<endl;

    //for(int i = 0;i < mAlignStart.size();i++){
    //		cout<<"\t"<<mReads[i].mName<<"\t"<<mAlignStart[i]<<endl;
    //	}
  }

  mNumBaseSeqs = mBaseSegments.size();

}


//----------------------------------------
// * WriteAce
//----------------------------------------
// Writes the contents of this contig out in 
// Ace format.  Some pieces of information are 
// not tracked and so not written out. 
// 
void 
Contig::WriteAce(ostream &fout)
{

  // ****  Sanity tests and fixups. *****

  // None of these conditions should happen, this is just a backup to ensure that
  // a consistent ace file is written out. 
  if ((mBaseSegments.size() == 0) || 
      (mNumBaseSeqs != mBaseSegments.size()) || 
      (mNumBaseSeqs ==0)){
    cerr<<"WARNING: "<<mName<<" mBaseSegments.size()="<<mBaseSegments.size();
    cerr<<" mNumBaseSeqs = "<<mNumBaseSeqs<<endl;
    ComputeSyntheticBS();
    cerr<<"\tFIXED: "<<mName<<" mBaseSegments.size()="<<mBaseSegments.size();
    cerr<<" mNumBaseSeqs = "<<mNumBaseSeqs<<endl;
  }

  if (mAlignStart.size() != mReads.size()){
    cerr<<"ERROR mAlignStart.size()="<<mAlignStart.size()<<" mReads.size()="<<mReads.size()<<endl;
  }

  // **** CO **************
  fout<<endl;
  
  // Ensure minimal sanity. 
  if (mPaddedBases != size()) mPaddedBases = size();
  if (mNumEntries != mReads.size()) mNumEntries=mReads.size();


  fout<<"CO "<<mName<<" "<<mPaddedBases<<" "<<mNumEntries<<" "<<mNumBaseSeqs<<" U"<<endl;
  for(int i = 0;i< (*this).size();i++){
    if ( (!(i%50)) && (i != 0)) fout<<endl;
    fout<<(*this)[i];
  }
  fout<<endl;
  fout<<endl;
  
  // **** BQ **************
  fout<<"BQ"<<endl;
  int breakCount = 0;

  for(int i = 0;i< mConsensusQuality.size();i++){

    // The Ace format only writes out quality values for the unpadded
    // bases. 
    if ((*this)[i] != '*'){
      breakCount++;
      fout<<mConsensusQuality[i]<<" ";
    }

    if ( (!(breakCount%50)) && (breakCount != 0)) {
      fout<<endl;
      breakCount++;
    }

  }
  fout<<endl<<endl;



  // ***** AF  **************
  for(int i = 0;i< mReads.size();i++){
    fout<<"AF "<<mReads[i].mName<<" "<<mReadSense[i]<<" "<<mAlignStart[i]<<endl;	 
  }	 


  // ***** BS  **************
  for(int i = 0;i< mBaseSegments.size();i++){
    fout<<"BS "<<mBaseSegments[i].mStart<<" "<<mBaseSegments[i].mEnd<<" "<<mBaseSegments[i].mName<<endl;
  }
  fout<<endl;


  // ***** RD/DS  **************

  for(int i = 0;i< mReads.size();i++){
    // KJD The 0 0 for whole_read_infos and read_tags is fudged.  
    fout<<"RD "<<mReads[i].mName<<" "<<mReads[i].size()<<" 0 0"<<endl;	 
    for(int j = 0;j < mReads[i].size();j++){
      if ( (!(j%50)) && (j != 0)) fout<<endl;
      fout<<mReads[i][j];
    }
    fout<<endl<<endl;

    // KJD QA Lines need to be computed for real... for now it's this horrible kludge. 
    if (mQALines.size() == 0){
      fout<<"QA "<<1<<" "<<mReads[i].size()-1<<" "<<1<<" "<<mReads[i].size()-1<<endl; // KJD Kludge default. 
    }else fout<<mQALines[i]<<endl; // KJD Unparsed but recorded. 

    fout<<"DS "<<mReads[i].mDescription<<endl;
    fout<<endl;
  }
  fout<<endl;
}
	 

//----------------------------------------
// * GetTrimmedDepadded
//----------------------------------------
//
void 
Contig::GetTrimmedDepadded(Sequence &Seq)
{
  Seq.mName = mName;
  Seq.mDescription = mDescription;
  for(int i = mTrimLeft;i<= mTrimRight;i++){
    if ((*this)[i] != '*'){
      Seq.push_back((*this)[i]);
    }
  }
}

//----------------------------------------
// * PrintBACTrimmedContig
//----------------------------------------
//
int
Contig::PrintBACTrimmedContig(ostream &fout,int minContigLength)
{
  int b;
  int pIdx = 0; // printed character index.  

  // See if contig meets minimum length and that there is some trim there.
  // could fail if the number of pads is very large.    
  if ((mBACTrimRight - mBACTrimLeft) >= minContigLength){
    if (mBACTrimLeft >=0){  

      fout<<">"<<mName<<" "<<mDescription<<endl;
		
      for(b = mBACTrimLeft;b<= mBACTrimRight;b++){
		
	// At the moment, pads are simply not allowed in the output.  All 
	// pads are simply skipped. 
	if ((*this)[b] != '*'){
	  fout<<(*this)[b];
	  if (((pIdx+1) % 50) == 0) fout<<endl;		  
	  pIdx++;
	}
      }
      // Except in the rare case where the line ended on an even 
      // multiple of 50, and an endl is already present, put in a 
      // line feed. 
      if ((pIdx % 50) != 0){
	fout<<"\n";
      }		
    }
  }
  return(pIdx);
}

//----------------------------------------
// * ChooseBestTemplatePair
//----------------------------------------
// Often there is more than one try at a template pair. 
// This results in more than two reads with the same template ID. 
// Sometimes there are four, representing a retry of the pair, 
// and sometimes there are three, representing a retry of one
// member of the pair.   This function attempts to pick the best
// pair out of whatever set of reads have the same template ID. 
//
void
Contig::ChooseBestTemplatePair(vector<int> &ReadIdxs,int &spanOut){
  
  vector<int> bugReads;
  bugReads = ReadIdxs;

  vector<int> UReads;
  vector<int> CReads;

  // Separate them into U reads and C reads
  for(int i = 0;i< ReadIdxs.size();i++){
    int r = ReadIdxs[i];
    int sense = mReadSense[r];
    //	 cerr<<mReads[r].mName<<"\t"<<mReadSense[r]<<endl;

    if (sense == 'U') UReads.push_back(r);
    else CReads.push_back(r);
  }

  ReadIdxs.clear();
  spanOut = 0;
  
  // If either one of these has size zero, then we don't have a 
  // pair.  
  if ((CReads.size() ==0) || (UReads.size() == 0)) return;

 
  // Now, for each U read, look for the best C read. 
  // 
  vector<int> bestMate,mateSpan;
  bestMate.resize(UReads.size());
  mateSpan.resize(UReads.size());
  for(int i = 0;i < mateSpan.size();i++) mateSpan[i] = 0;


  for(int u = 0;u < UReads.size();u++){
    int uIdx = UReads[u];
	 
    int maxSpan = -9999,maxReadIdx=0;
    for(int c = 0;c< CReads.size();c++){
      int cIdx = CReads[c];
		
      int spanStart = mAlignStart[uIdx] + mReads[uIdx].size();
      int spanEnd = mAlignStart[cIdx]-1;
      int span = spanEnd - spanStart;

      if (span > maxSpan){
	maxSpan = span;
	maxReadIdx = cIdx;
      }
    }

    bestMate[u] = maxReadIdx;
    mateSpan[u] = maxSpan;
  }  

  // Now find the best overall -pair-
  int maxSpan = -9999,bestIdx = -1;
  for(int u = 0;u < UReads.size();u++){
    if (mateSpan[u] > maxSpan) {
      maxSpan = mateSpan[u];
      bestIdx = u;
    }
  }

  if ((bestIdx < 0) || (bestIdx > UReads.size())){
    return;
  }

  // Return ReadIdxs with only the best pair. 
  ReadIdxs.push_back(UReads[bestIdx]);
  ReadIdxs.push_back(bestMate[bestIdx]);

  spanOut = maxSpan;


  // KJD BUG figure out what is happening
  // KJD     when have negative span==> appears that there are just a 
  // lot of pairs with no gap between them. 
  /*
    if (maxSpan < 0){
    for(int i = 0;i< bugReads.size();i++){
    int r = bugReads[i];

    cout<<r<<"\t";
    cout<<mReads[r].mName<<"\t";
    cout<<mReadSense[r]<<"\t";
    cout<<mAlignStart[r]<<"\t";
    cout<<mReads[r].size()<<"\t";
    if (mReadSense[r] == 'U') cout<<"start="<<mAlignStart[r]+mReads[r].size()<<endl;
    else cout<<"end="<<mAlignStart[r]<<endl;
    }
    }
  */
}




//----------------------------------------
// * GetTemplateInfo
//----------------------------------------
// Gathers information about all of the templates in the contig, computing 
// template size and other stats for each template pair. Basically repackages
// the information that is already available in an easy to use data structure. 
// 
void
Contig::GetTemplateInfo(vector<TemplateInfo> &TInfo)
{
  // KJD make sure there isn't already some template info in the list...
  TInfo.clear();

  int bugSize = 0;
  if (mTemplate2ReadMap.size() == 0) return;
  else bugSize = mTemplate2ReadMap.size();


  vector<string> keys;
  mTemplate2ReadMap.GetKeys(keys);

  // for each template id
  for(int i = 0;i< keys.size();i++){
    int span = -9999;
		
    TemplateInfo temp;
    temp.mTID = keys[i];
		
    vector<int> Reads;
    mTemplate2ReadMap.GetValues(temp.mTID,Reads);
		
    // Get the reads that go with that template ID
    // Often there are more than two reads with a given template ID, so we have
    // to decide which pair to use in that case. 
    if (Reads.size() > 2){
      ChooseBestTemplatePair(Reads,span);
			
      // If this fails, then skip this template.  KJD Why would this ever happen?
      if (Reads.size() != 2) {
	cout<<"ERROR: After ChooseBestTemplatePair, there are not 2 reads."<<endl;
	continue;			
      }
    }
	 
    if (Reads.size() == 0){
      cerr<<" ERROR: No reads for template."<<endl;
      continue;
    }

    string desc = mReads[Reads[0]].mDescription;
    int idx = desc.find("INSERT_SIZE:");
    temp.mInsertSize = 0;
    if (idx >=0){
      idx += sizeof("INSERT_SIZE:"); // point to space after token. 
      temp.mInsertSize = NextInt(desc,idx);
    }

    // Assume at most two for size computation. 
    int firstReadIdx = Reads[0];
    int secondReadIdx = -1;

    // If we have a clean pair. 
    if (Reads.size() == 2){		
      temp.mbPair = true;

      // Figure out which read is earlier in the contig. 
      if (mAlignStart[Reads[0]] <= mAlignStart[Reads[1]]){
	firstReadIdx = Reads[0];
	secondReadIdx = Reads[1];
      }else{
	firstReadIdx = Reads[1];
	secondReadIdx = Reads[0];
      }
		
      // Since there is a clean pair we can compute the size explicitly. 
      temp.mLength = mAlignStart[secondReadIdx] - mAlignStart[firstReadIdx]+mReads[secondReadIdx].size();
      temp.mRSize1 = mReads[firstReadIdx].size();
      temp.mRSize2 = mReads[secondReadIdx].size();
    }else{
      // Not a pair, so the length is simply the size of the single read.  
      temp.mbPair = false;
      firstReadIdx = Reads[0];
      temp.mLength = mReads[firstReadIdx].size(); 
      temp.mRSize1 = mReads[firstReadIdx].size(); 
    }
	
    // A little sanity check. 
    //if (span != -9999){
    //	if (span != temp.mLength) {
    //	  cerr<<"BUG: span="<<span<<"\tmLength="<<temp.mLength<<endl;
    //	}
    // }

    // read info for each read. 
    if (mReadTemplateDir[firstReadIdx] == 0) temp.mTDir1 ="fwd";
    else temp.mTDir1="rev";
	 
    // KJD Better would be just to save the read Idx in the template info. 
    //cout<<"mReads.size()="<<mReads.size()<<" firstReadIdx="<<firstReadIdx<<" secondReadIdx="<<secondReadIdx<<endl;

    temp.mRName1 = mReads[firstReadIdx].mName;
    temp.mRDir1 = mReadSense[firstReadIdx];
    temp.mRStart1 = mAlignStart[firstReadIdx];
	 
    if (Reads.size() > 1){		
		
      if (mReadTemplateDir[secondReadIdx] == 0) temp.mTDir2 ="fwd";
      else temp.mTDir2="rev";
		
      temp.mRName2 = mReads[secondReadIdx].mName;		
      temp.mRDir2 = mReadSense[secondReadIdx];
      temp.mRStart2 = mAlignStart[secondReadIdx];		

      // Now check if pair is in the right orientation. To be in 
      // the right orientation, the downstream read must be on 
      // the uncomplemented strand (relative to the assembled 
      // contig) and the upstream read must be on the complemented
      // strand.  It doesn't matter whether the read is fwd/rev end. 
      if ((temp.mRDir1 == 'U') && (temp.mRDir2 == 'C')){
	temp.mbOrientationOK = true;
      }else{
	temp.mbOrientationOK = false;
      }
    }
    TInfo.push_back(temp);
  }
}

//----------------------------------------
// * BuildBadTemplateStartEnd
//----------------------------------------
// Given the TemplateInfo for the contig, this function builds a 
// 
// 
// KJD NOTE: starts that occur before the beginning of the 
// contig are all mapped to contig position zero. 
//
void
Contig::BuildBadTemplateStartEnd(vector<TemplateInfo> &TInfo,
				 Int2IntMap &startmap,Int2IntMap &endmap)
{ 
  for(int i = 0;i< TInfo.size();i++){


    // If it's a part of a pair, and the orientation is 
    // incorrect, record the start and end of the template 
    // covered region. 
    if (TInfo[i].mbPair && !TInfo[i].mbOrientationOK){
      int startCIdx = TInfo[i].mRStart1;

      // Map negative contig indices to zero. 
      if (startCIdx < 0) startCIdx = 0;

      if (startmap.Contains(startCIdx)){
	startmap[startCIdx] = startmap[startCIdx]+1;
      }else{
	startmap[startCIdx] = 1;
      }		  
		
      int endCIdx = TInfo[i].mRStart2 + TInfo[i].mRSize2; // KJD 
      if (endCIdx < 0) endCIdx = 0;

      if (endmap.Contains(endCIdx)){
	endmap[endCIdx] = endmap[endCIdx]+1;
      }else{
	endmap[endCIdx] = 1;
      }
    }		
  }
}

//----------------------------------------
// * BuildTemplateStartEnd
//----------------------------------------
// Given the TemplateInfo for the contig, this function builds a 
// 
// 
// KJD NOTE: starts that occur before the beginning of the 
// contig are all mapped to contig position zero. 
//
void
Contig::BuildTemplateStartEnd(vector<TemplateInfo> &TInfo,
			      Int2IntMap &startmap,Int2IntMap &endmap)
{ 
  for(int i = 0;i< TInfo.size();i++){


    // If it's a part of a pair, and the orientation is 
    // correct, record the start and end of the template 
    // covered region. 
    if (TInfo[i].mbPair && TInfo[i].mbOrientationOK){
      int startCIdx = TInfo[i].mRStart1;

      // Map negative contig indices to zero. 
      if (startCIdx < 0) startCIdx = 0;

      if (startmap.Contains(startCIdx)){
	startmap[startCIdx] = startmap[startCIdx]+1;
      }else{
	startmap[startCIdx] = 1;
      }		  
		
      int endCIdx = TInfo[i].mRStart2 + TInfo[i].mRSize2; // KJD 
      if (endCIdx < 0) endCIdx = 0;

      if (endmap.Contains(endCIdx)){
	endmap[endCIdx] = endmap[endCIdx]+1;
      }else{
	endmap[endCIdx] = 1;
      }
    }		
  }
}


//----------------------------------------
// * BuildTemplateSpanStartEnd
//----------------------------------------
// For this function, start is one after the last base of the left read
// end is one before the first base of the right read in a pair. 
//
void
Contig::BuildTemplateSpanStartEnd(vector<TemplateInfo> &TInfo,
				  Int2IntMap &startmap,Int2IntMap &endmap)
{ 
  for(int i = 0;i< TInfo.size();i++){
	 
	 
    // If it's a part of a pair, and the orientation is 
    // correct, record the start and end of the template 
    // covered region. 
    if (TInfo[i].mbPair && TInfo[i].mbOrientationOK){
      int startCIdx = TInfo[i].mRStart1+TInfo[i].mRSize1;

      // Map negative contig indices to zero. 
      if (startCIdx < 0) startCIdx = 0;

      if (startmap.Contains(startCIdx)){
	startmap[startCIdx] = startmap[startCIdx]+1;
      }else{
	startmap[startCIdx] = 1;
      }		  
		
      //		int endCIdx = TInfo[i].mRStart2 + TInfo[i].mRSize2; // KJD 

      int endCIdx = TInfo[i].mRStart2 -1;

      if (endCIdx < 0) endCIdx = 0;

      if (endmap.Contains(endCIdx)){
	endmap[endCIdx] = endmap[endCIdx]+1;
      }else{
	endmap[endCIdx] = 1;
      }
    }		
  }
}

//----------------------------------------
// * GetTemplateSpan
//----------------------------------------
// Sort of like the template coverage, but counts the number
// of templates which -span- any given base.  
//
void
Contig::GetTemplateSpan(vector<int> &TCover)
{
  Int2IntMap starts,ends;

  TCover.resize(size());
  for(int i = 0;i< TCover.size();i++) TCover[i]=0;

  mTInfo.erase(mTInfo.begin(),mTInfo.end());
  GetTemplateInfo(mTInfo);


  BuildTemplateSpanStartEnd(mTInfo,starts,ends);

  int coverage = 0;
  for(int cIdx = 0;cIdx < size();cIdx++){
    // Add any new coverage that appears at this 
    // index. 
    if (starts.Contains(cIdx)){
      coverage = coverage + starts[cIdx];
      //cout<<"starts["<<cIdx<<"]="<<starts[cIdx]<<endl;
    }
    // Record the coverage at this location. 
    TCover[cIdx] = coverage;

    // Subtract coverage of any template pairs that end
    // at this location. 
    if (ends.Contains(cIdx)){
      coverage = coverage - ends[cIdx];
      //cout<<"ends["<<cIdx<<"]="<<ends[cIdx]<<endl;
    }
  }
}

//----------------------------------------
// * GetBadTemplateCoverage
//----------------------------------------
// Returns the bad template coverage at each contig position. 
// This is only crudely done. 
//
void
Contig::GetBadTemplateCoverage(vector<int> &TCover)
{
  Int2IntMap starts,ends;

  TCover.resize(size());
  for(int i = 0;i< TCover.size();i++) TCover[i]=0;

  mTInfo.erase(mTInfo.begin(),mTInfo.end());
  GetTemplateInfo(mTInfo);

  
  BuildBadTemplateStartEnd(mTInfo,starts,ends);
  
  int coverage = 0;
  for(int cIdx = 0;cIdx < size();cIdx++){
    // Add any new coverage that appears at this 
    // index. 
    if (starts.Contains(cIdx)){
      coverage = coverage + starts[cIdx];
      //cout<<"starts["<<cIdx<<"]="<<starts[cIdx]<<endl;
    }
    // Record the coverage at this location. 
    TCover[cIdx] = coverage;

    // Subtract coverage of any template pairs that end
    // at this location. 
    if (ends.Contains(cIdx)){
      coverage = coverage - ends[cIdx];
      //cout<<"ends["<<cIdx<<"]="<<ends[cIdx]<<endl;
    }
  }
}

//----------------------------------------
// * GetTemplateCoverage
//----------------------------------------
// Returns the template coverage at each contig position. 
//
void
Contig::GetTemplateCoverage(vector<int> &TCover)
{
  Int2IntMap starts,ends;
	
  int cSize = (*this).size();
  //	cout<<"cSize="<<cSize<<endl;
  if (cSize == 0) return;

  TCover.resize(cSize);
  for(int i = 0;i< TCover.size();i++) TCover[i]=0;

  vector<TemplateInfo> TInfo;
  GetTemplateInfo(TInfo);

  BuildTemplateStartEnd(TInfo,starts,ends);
  
  int coverage = 0;
  for(int cIdx = 0;cIdx < size();cIdx++){
    // Add any new coverage that appears at this 
    // index. 
    if (starts.Contains(cIdx)){
      coverage = coverage + starts[cIdx];
      //cout<<"starts["<<cIdx<<"]="<<starts[cIdx]<<endl;
    }
    // Record the coverage at this location. 
    TCover[cIdx] = coverage;

    // Subtract coverage of any template pairs that end
    // at this location. 
    if (ends.Contains(cIdx)){
      coverage = coverage - ends[cIdx];
      //cout<<"ends["<<cIdx<<"]="<<ends[cIdx]<<endl;
    }
  }
}

//----------------------------------------
// * CalcP20Coverage
//----------------------------------------
//
int 
Contig::CalcP20Coverage(int idx)
{
  int coverage = 0;
  for(int r = 0;r< mReads.size();r++){
    int readIdx = CIdx2RIdx(idx,r);
    if ((readIdx>=0) && (readIdx<mReadsQuality[r].size())){
      if (mReadsQuality[r][readIdx] >= 20){
	coverage++;
      }
    }
  }
  return(coverage);
}


//----------------------------------------
// * CalcCoverage
//----------------------------------------
//
int 
Contig::CalcCoverage(int idx)
{
  int coverage = 0;
  for(int r = 0;r< mReads.size();r++){
    int readIdx = CIdx2RIdx(idx,r);
    if ((readIdx>=0) && (readIdx<mReads[r].size())){
      coverage++;
    }
  }
  return(coverage);
}


//----------------------------------------
// * GetCoverage
//----------------------------------------
//
void 
Contig::GetCoverage(vector<int> &Coverage)
{
  Coverage.resize(size());
  
  for(int i = 0;i< size();i++){
    Coverage[i] = CalcCoverage(i);
  }
}


//----------------------------------------
// * GetHQDescrepency
//----------------------------------------
//
void 
Contig::GetHQDescrepency(vector<int> &HQDesc)
{

  //cout<<"GetHQ descrepency. size="<<size()<<endl;

  // OK, right now it is just descrepency. 
  if (mReadsQuality.size() != mReads.size()){
    cerr<<"ReadQuality doesn't match number of reads."<<endl;
    return;
  }

  HQDesc.resize(size());
	 
  for(int i = 0;i< HQDesc.size();i++) HQDesc[i] = 0;
	 
  for(int r = 0;r< mReads.size();r++){
    for(int rIdx = 0;rIdx < mReads[r].size();rIdx++){
      int cIdx = RIdx2CIdx(rIdx,r);
				
      // If it's a valid contig index....
      if ((cIdx >=0) && (cIdx < size())){		
	// If base does not match consensus 
	char base = toupper(mReads[r][rIdx]);
	char consensus = toupper((*this)[cIdx]);
	int qual = mReadsQuality[r][rIdx];
	if (qual >= 20){
	  if (base != consensus) {
	    HQDesc[cIdx]++;
	  }
	}
      }else{
	// cerr<<"Warining cIdx="<<cIdx<<" size="<<size()<<endl;
      }
    }
  }
}


//----------------------------------------
// * GetWGSBACCoverage
//----------------------------------------
//
void
Contig::GetWGSBACCoverage(vector<int> &WGSCover,vector<int> &BACCover)
{
  int cSize = (*this).size();
  WGSCover.resize(cSize);
  BACCover.resize(cSize);

  for(int i = 0;i< cSize;i++){
    WGSCover[i] = 0;
    BACCover[i] = 0;
  }

  bool bWGS;

  for(int r = 0;r < mReads.size();r++){
    if (Contains(mReads[r].mName,"gnl")) bWGS = true;
    else bWGS = false;

    int rSize = mReads[r].size();
    for(int rIdx = 0;rIdx < rSize;rIdx++){
      int cIdx = RIdx2CIdx(rIdx,r);
      if ((cIdx >=0) && (cIdx < cSize)){
	if (bWGS) WGSCover[cIdx]++;
	else BACCover[cIdx]++;
      }
    }
  }
}


//----------------------------------------
// * PrintTemplateInfo
//----------------------------------------
//
void
Contig::PrintTemplateInfo(ostream &fout)
{

  vector<TemplateInfo> TInfo;

  GetTemplateInfo(TInfo);

  double totalInsertSize = 0;

  // for each template id
  for(int i = 0;i< TInfo.size();i++){
    fout<<TInfo[i].mTID;	 
    fout<<" ("<<TInfo[i].mLength<<")\t";
    totalInsertSize+=TInfo[i].mLength;

    fout<<TInfo[i].mRName1<<"("<<TInfo[i].mTDir1<<","<<TInfo[i].mRDir1<<")\t";
    fout<<TInfo[i].mRStart1<<"\t"<<TInfo[i].mRSize1<<"\t";

    if (TInfo[i].mRName2 != ""){
      fout<<TInfo[i].mRName2<<"("<<TInfo[i].mTDir2<<","<<TInfo[i].mRDir2<<")\t";
      fout<<TInfo[i].mRStart2<<"\t"<<TInfo[i].mRSize2<<"\t";
    }
    fout<<endl;
  }
  double avgInsertSize = (double)totalInsertSize/(double)TInfo.size();
  fout<<endl;
  fout<<"Total pairs: "<<TInfo.size()<<"\tAverage insert size: "<<avgInsertSize<<endl;
  fout<<endl;
}

//----------------------------------------
// * RecordTemplateInfo
//----------------------------------------
// Extract the template2ReadMap, Read2TemplateMap,
// and the ReadTemplateDir from the description lines of
// contigs.   
// 
// KJD Should probably just create a TInfo vec from the 
// KJD start.  That is, bundle GetTemplateInfo with this. 
// 
// 
void 
Contig::RecordTemplateInfo(int ReadIdx)
{
  string desc = mReads[ReadIdx].mDescription;
  int dirInt = 0;
  string templateID;

  int idx = desc.find("TEMPLATE:");
  if (idx >=0){
    idx += sizeof("TEMPLATE:"); // point to space after token. 
    templateID = GetNextWSDelimitedString(desc,idx,idx);

    if (!ContainsOnlyWS(templateID)){
      idx = desc.find("DIRECTION:");
      if (idx >=0){
	idx+= sizeof("DIRECTION:");
	string dir = GetNextWSDelimitedString(desc,idx,idx);	
	if (Contains(dir,"fwd")) dirInt = 0;
	else dirInt = 1;
      }
      mTemplate2ReadMap.Add(templateID,ReadIdx);
      mRead2TemplateMap.Add(mReads[ReadIdx].mName,templateID);
      mReadTemplateDir.Add(ReadIdx,dirInt);  
    }
  }
}


//----------------------------------------
// * CalcAverageTrimmedCoverage
//----------------------------------------
// Compute the average coverage for this contig. 
// Coverage is computed over the trimmed portion of 
// reads.   Requires CalculateReadTrims to be called
// first. 
// 
float
Contig::CalcAverageTrimmedCoverage()
{
  long totalC=0;
  for(int c = 0;c< size();c++){
    int bcover = CalcCoverageTrimmed(c);
    totalC+=bcover;
  }

  float retval = (float)totalC/(float)size();
  return(retval);
}


//----------------------------------------
// * CalcAverageCoverage
//----------------------------------------
// Compute the average coverage for this contig. 
// Coverage is computed over entire read. 
// 
float
Contig::CalcAverageCoverage()
{
  long totalC=0;
  for(int c = 0;c< size();c++){
    int bcover = CalcCoverage(c);
    totalC+=bcover;
  }

  float retval = (float)totalC/(float)size();
  return(retval);
}


long 
Contig::CountP20Bases()
{
  long total = 0;
  for(int c = 0;c < size();c++){
    int bcover = CalcP20Coverage(c);
    total+=bcover;
  }
  return(total);
}


//----------------------------------------
// * CalcAverageP20Coverage
//----------------------------------------
// Compute the average coverage for this contig. 
// Coverage is computed over entire read. 
// 
float
Contig::CalcAverageP20Coverage()
{
  long totalC=0;
  for(int c = 0;c< size();c++){
    int bcover = CalcP20Coverage(c);
    totalC+=bcover;
  }

  float retval = (float)totalC/(float)size();
  return(retval);
}



//----------------------------------------
// * CalcCoverageTrim
//----------------------------------------
//
int 
Contig::CalcCoverageTrimmed(int idx)
{
  int coverage = 0;
  for(int r = 0;r< mReads.size();r++){	 
    int readIdx = CIdx2RIdx(idx,r);	 

    if ((readIdx >= 0) && 
	(readIdx < mReads[r].size()) &&
	(readIdx >= mReadTrimLeft[r])&&
	(readIdx <= mReadTrimRight[r])){
      // And, finally, ignore it if it's marked with soft trimming. 
      //if (isupper(mReads[r][readIdx])){
      //  coverage++;
      //	}
      coverage++;
    }
  }
  return(coverage);
}

//----------------------------------------
// * GetMatchingBasesTrimmed
//----------------------------------------
// Figures out how many bases of the specified 
// read match bases of the contig consensus. 
// Only examines bases inside the contig and 
// read trim regions.  
// 
int
Contig::GetMatchingBasesTrimmed(int readIdx)
{
  int matchCount = 0;

  // For each trimmed base of the contig
  for(int cIdx = mTrimLeft;cIdx <= mTrimRight;cIdx++){
    // Find the corresponding base in the read...
    int readBaseIdx = CIdx2RIdx(cIdx,readIdx);

    // Skip bases outside the read or read trim
    if ((readBaseIdx >= 0) && 
	(readBaseIdx < mReads[readIdx].size()) &&
	(readBaseIdx >= mReadTrimLeft[readIdx])&&
	(readBaseIdx <= mReadTrimRight[readIdx])){

      char contigBase = tolower((*this)[cIdx]);
      char readBase = tolower(mReads[readIdx][readBaseIdx]);

      // skip pads in the contig...
      if (contigBase != '*'){
	if (contigBase == readBase){
	  matchCount++;
	}
      }
    }
  }

  return(matchCount);
}

//----------------------------------------
// * GetMatchingBases
//----------------------------------------
// Figures out how many bases of the specified 
// read match bases of the contig consensus. 
// Only examines bases inside the contig trim. 
// 
int
Contig::GetMatchingBases(int readIdx)
{
  int matchCount = 0;

  // For each base in the contig.  
  for(int cIdx = 0;cIdx <= size();cIdx++){
    // Find the corresponding base in the read...
    int readBaseIdx = CIdx2RIdx(cIdx,readIdx);

    // Skip bases outside the read 
    if ((readBaseIdx >= 0) && 
	(readBaseIdx < mReads[readIdx].size())){
      char contigBase = tolower((*this)[cIdx]);
      char readBase = tolower(mReads[readIdx][readBaseIdx]);
		
      // skip pads in the contig...
      if (contigBase != '*'){
	if (contigBase == readBase){
	  matchCount++;
	}
      }
    }
  }

  return(matchCount);
}


//----------------------------------------
// * HasMateInContig
//----------------------------------------
// Tests whether the read specified by read name has a 
// mate in this same contig. 
// 
bool
Contig::HasMateInContig(string R1)
{
  if (mRead2TemplateMap.Contains(R1)){
    string TemplateID = mRead2TemplateMap[R1];
	 
    vector<int> ReadIdxs;
    mTemplate2ReadMap.GetValues(TemplateID,ReadIdxs);
    if (ReadIdxs.size() > 1) return(true);
    else return(false);
  }else{
    cerr<<R1<<" does not have a template recorded in this contig."<<endl;
    return(false);
  }
}


//----------------------------------------
// * HasMateInContig
//----------------------------------------
// Tests whether the read specified by read name has a 
// mate in this same contig. 
// 
bool
Contig::HasGoodMateInContig(string R1)
{
  if (mRead2TemplateMap.Contains(R1)){
    string TemplateID = mRead2TemplateMap[R1];
	 
    vector<int> ReadIdxs;
    mTemplate2ReadMap.GetValues(TemplateID,ReadIdxs);
	 
    for(int i = 0;i<mTInfo.size();i++){
      if (mTInfo[i].mTID == TemplateID){
	if (mTInfo[i].mbPair){
	  if (mTInfo[i].mbOrientationOK) return(true);
	  else return(false);
	}else return(false);
      }
    }
  }
  return(false);
}


//  mTemplate2ReadMap.Add(templateID,ReadIdx);
//  mRead2TemplateMap.Add(mReads[ReadIdx].mName,templateID);
//  mReadTemplateDir.Add(ReadIdx,dirInt);  
//  Contig::RecordTemplateInfo(int ReadIdx)

//----------------------------------------
// * RemoveRead
//----------------------------------------
// Removes all traces of a read from the contig. 
// 
void
Contig::RemoveRead(string readName)
{

  if (!mReadIdxMap.Contains(readName)){
    cerr<<"Can't remove "<<readName<<" Not in index."<<endl;
  }

  int rIdx = mReadIdxMap[readName];

  // Sanity check. 
  if (rIdx >= mReads.size()){
    cerr<<"ERROR: Remove rIdx ="<<rIdx<<" mReads.size="<<mReads.size()<<endl;
  }

  // Bit of a mess to remove from this so skip it for now. 
  // mTemplate2ReadMap;
  
  mReads.erase(mReads.begin()+rIdx);
  mReadsQuality.erase(mReadsQuality.begin()+rIdx);
  mAlignStart.erase(mAlignStart.begin()+rIdx);

  mReadSense.erase(mReadSense.begin()+rIdx);
  mReadTrimLeft.erase(mReadTrimLeft.begin()+rIdx);
  mReadTrimRight.erase(mReadTrimRight.begin()+rIdx);

  // Since it is a mess to modify the base segments, I'll 
  // just erase them all for now.  
  mBaseSegments.erase(mBaseSegments.begin(),mBaseSegments.end());
  // Add in one dummy base segment, just to not break some parsers
  /*
    mBaseSegments.resize(1);
    mBaseSegments[0].mStart=0;
    mBaseSegments[0].mEnd=0;
    mBaseSegments[0].mName = "NotPhrapConsensus";
  */

	
  mQALines.erase(mQALines.begin()+rIdx);
  
  // Rebuild the read index map 

  mReadIdxMap.erase(mReadIdxMap.begin(),mReadIdxMap.end());
  for(int r = 0;r< mReads.size();r++){
    mReadIdxMap[mReads[r].mName] = r;
  }

  // Erase all of the template info and rebuild it.  
  mTemplate2ReadMap.clear();
  mRead2TemplateMap.clear();
  mReadTemplateDir.clear();
  for(int r = 0;r<mReads.size();r++){
    RecordTemplateInfo(r);
  }
  mNumEntries--;
}


//----------------------------------------
// * RebuildTemplateInfo
//----------------------------------------
//
void 
Contig::RebuildTemplateInfo()
{

  // Rebuild the read index map 

  mReadIdxMap.erase(mReadIdxMap.begin(),mReadIdxMap.end());
  for(int r = 0;r< mReads.size();r++){
    mReadIdxMap[mReads[r].mName] = r;
  }

  // Erase all of the template info and rebuild it.  
  mTemplate2ReadMap.clear();
  mRead2TemplateMap.clear();
  mReadTemplateDir.clear();
  for(int r = 0;r<mReads.size();r++){
    RecordTemplateInfo(r);
  }
}


//----------------------------------------
// * RemoveTails
//----------------------------------------
// Removes "unanchored" reads in the tails of contigs specified 
// by left/right contig indices.   An anchored read is a read 
// that has a mate within the contig.  An unanchored read does
// not have a mate within the contig.   
//
void 
Contig::RemoveTails(int leftTail,int rightTail,vector<string> &RemovedReads)
{
  vector<string> KillSet;

  // Find the reads 
  StringSet OffendingReads;
  int readsSize = mReads.size();
  for(int r = 0;r < mReads.size();r++){
    // Compute the effective read index based on the contig index. 
    int lReadIdx = CIdx2RIdx(leftTail,r);
    int rReadIdx = CIdx2RIdx(rightTail,r);

    // If the read index of the left tail point is greater than zero, 
    // then the left tail falls inside, or to the right of the read, making
    // this an offending read. 
    if (lReadIdx > 0){
      //cout<<"\tLOffender"<<flush;
      OffendingReads.Add(mReads[r].mName);
      if (!HasGoodMateInContig(mReads[r].mName)){
	//cout<<"\t NoMate"<<flush;
	KillSet.push_back(mReads[r].mName);
      }else{
	//cout<<"\t HasMate"<<flush;
      }
    }

    int a = rReadIdx;
    int b = mReads[r].size();
    //cout<<"rrIdx="<<a<<" size="<<b<<flush;
    if (a < b){
      //cout<<" rrIdx<size "<<flush;
      //cout<<"\tROffender"<<flush;
      OffendingReads.Add(mReads[r].mName);
      if (!HasGoodMateInContig(mReads[r].mName)){
	KillSet.push_back(mReads[r].mName);
	//cout<<"\t NoMate"<<flush;
      }else{
	//cout<<"\t HasMate"<<flush;
      }
    }
    // cout<<endl;
  }

  //  cout<<endl<<"Total Offending Reads: "<<OffendingReads.size()<<flush;

  //cout<<"\tRemoveReads: "<<KillSet.size()<<endl;

  
  for(int i = 0;i< KillSet.size();i++){	 
    RemoveRead(KillSet[i]);
    RemovedReads.push_back(KillSet[i]); // Make global list of killed reads. 
  }

  // Recompute the consensus in the region of interest. 

  ComputeConsensus(rightTail,size()-1);

  // Actually remove the pad tail. 
  TrimPadTail();
}


//----------------------------------------
// * FindReadSets
//----------------------------------------
// Given two contig positions, it returns the set of 
// reads which are between those positions and unanchored 
// to either side by a template pair (MiddleSet), those
// reads which are to the left of the contig position or 
// anchored to the left by a template pair (LeftSet), 
// and the set of reads which are to the right of the contig
// position or anchored to the right by a template pair (RightSet). 
//
// This function was written to aid in splitting contigs 
// at template crashes.  
// 

// There is no longer just a leftmost and rightmost altered, 
// but four points.  Seems confusing now to have that in this
// function. 
//
//
// KJD 
// KJD This must handle conditions where crash goes to left 
// KJD or right end or covers entire contig. 
// KJD 
//
void 
Contig::FindReadSets(int leftCrash, int rightCrash, 
		     vector<string> &LeftVec,
		     vector<string> &MiddleVec, 
		     vector<string> &RightVec)
{
  LeftVec.clear();
  MiddleVec.clear();
  RightVec.clear();

  // For a little sanity checking. 
  StringSet LeftSet,MiddleSet,RightSet;

  // Find the reads that fall in leftCrash/rightCrash interval.  
  StringSet OffendingReads;
  int readsSize = mReads.size();
  for(int r = 0;r < mReads.size();r++){
    // Convert read start/end into contig coordinates.  
    string RName = mReads[r].mName;
    int rSize = mReads[r].size();
    int rEnd = RIdx2CIdx(rSize-1,r);
    int rStart = RIdx2CIdx(0,r);  // KJD Important safety tip: rStart can be negative!!!
	 
    // KJD: Strictly speaking, this should be HasGoodMateInRight region
    // KJD: However, since we know that leftCrash/rightCrash mark a region 
    // KJD: of zero template coverage, we know the mate can't be from the opposite
    // KJD: side, if it exists. 
    bool bHasGoodMate = HasGoodMateInContig(RName);

    bool bFirmlyLeft = (leftCrash > rEnd);
    bool bOverlapsLeftCrash = ((rStart <= leftCrash) && (leftCrash <= rEnd));
    bool bUnanchoredLeft = (bOverlapsLeftCrash && (!bHasGoodMate));
    bool bAnchoredLeft = (bOverlapsLeftCrash && bHasGoodMate);
	 
    bool bFirmlyRight = (rightCrash < rStart);
    bool bOverlapsRightCrash = ((rStart <= rightCrash) && (rightCrash <= rEnd));
    bool bUnanchoredRight = (bOverlapsRightCrash && (!bHasGoodMate));
    bool bAnchoredRight = (bOverlapsRightCrash && (bHasGoodMate));
	 
    bool bBetweenLeftAndRightCrash = ((leftCrash < rStart) && (rEnd < rightCrash));
	 
    if (bUnanchoredLeft || bBetweenLeftAndRightCrash || bUnanchoredRight){
      MiddleVec.push_back(RName);
      //MiddleSet.Add(RName);
    }else{

      // LLLLLLLLLMMMMMMMMMMMMRRRRRRRRRRR
      //         XXXXXXXXXXXXXXXXX
      //
      // A read like XXXX, in order to go to the right side would have
      // to have a mate pair in the contig.  If the mate pair was on the 
      // left, MMMM wouldn't be a template crash.  Ditto if it is on the 
      // right.   Therefore, the existence of this condition would be a
      // bug. 
      // 
      if (bFirmlyRight || bAnchoredRight){
	RightVec.push_back(RName);
	//RightSet.Add(RName);
      }else{
	if (bFirmlyLeft || bAnchoredLeft){
	  LeftVec.push_back(RName);
	  //LeftSet.Add(RName);
	}else{
	  cerr<<"\tERROR BUG  This read has no home: "<<RName<<endl;
	}
      }
    }
  }
  
  // Sanity:  Every read should be accounted for. 
  int totalSize = LeftVec.size()+MiddleVec.size()+RightVec.size();
  if (totalSize != mReads.size()){
    cerr<<"ERROR: Left+Middle+Right="<<totalSize<<" mReads.size="<<mReads.size()<<endl;
    cerr<<"ERROR: Left = "<<LeftVec.size()<<"\tMiddle="<<MiddleVec.size()<<"\tRight="<<RightVec.size()<<endl;
  }
}

//----------------------------------------
// * CopyConsensus
//----------------------------------------
// Copies the consensus and consensus quality from one contig 
// to this one.  This contig is resized to accomodate the copied 
// sequence.  The copy starts at targetStart in this consensus, 
// and ranges from sourceStart to sourceEnd in the source contig. 
// 
void
Contig::CopyConsensus(int targetStart,Contig &C,int sourceStart,int sourceEnd)
{
  int copySize = sourceEnd - sourceStart +1;
  int minContigSize = copySize + targetStart;

  resize(minContigSize);
  mConsensusQuality.resize(minContigSize);

  int tIdx = targetStart;
  for(int cIdx = sourceStart;cIdx <=sourceEnd;cIdx++){
    (*this)[tIdx] = C[cIdx];
    mConsensusQuality[tIdx] = C.mConsensusQuality[cIdx];
    tIdx++;
  }
}


//----------------------------------------
// * FixupAlignStarts
//----------------------------------------
// A bit of a hack. When reads are copied to a new contig, 
// their alignStarts are relative to their source. 
// This renormalizes them relative to newStart, and returns
// the leftmost and rightmost extent of the reads in the 
// old coordinates.  
//
void 
Contig::FixupAlignStarts(int newStart,int &oldLeftMost,int &oldRightMost){

  // All of the alignStart values for this are with respect to the
  // source contig.   We're essentially starting from scratch, so
  // first determine the leftmost and rightmost bounds in the old 
  // coordinate system....
  oldLeftMost = 99999999;
  oldRightMost = 0;
  for(int r = 0;r < mReads.size();r++){
    if (mAlignStart[r] < oldLeftMost) oldLeftMost = mAlignStart[r];
    if (mAlignStart[r] > oldRightMost ) oldRightMost = mAlignStart[r];
  }

  // KJD A bit strained to store this, but we need it later and 
  // we're several functions deep, so hey...
  mOldLeftMost = oldLeftMost;
  mOldRightMost = oldRightMost;

  // OK, leftmostAlignStart is now the newStart base position, 
  // so renumber all of the alignStart locations
  for(int r = 0;r< mReads.size();r++){
    mAlignStart[r] = mAlignStart[r] - oldLeftMost + newStart;
  }
}

//----------------------------------------
// * TrimPadTail
//----------------------------------------
// Removing a bunch of reads in a split leaves a long 
// tail of pads.   This method trims that tail off. 
// 
void
Contig::TrimPadTail()
{
  //  cout<<"TrimPadTail. Contig.size="<<size()<<endl;

  Contig TempC;

  // If, for some crazy reason, the contig size is zero
  // there is nothing to do. 
  if (size() <= 0) return;
  
  int cIdx = size()-1;
  int numStar = 0;

  if ((*this)[cIdx] != '*') return;

  while((*this)[cIdx] == '*'){
    numStar++;
    cIdx--;
    if (cIdx < 0) break;
  }

  int newSize = size()-numStar;

  TempC = (*this);

  clear();
  mConsensusQuality.clear();
  //cout<<" C cleared size="<<size()<<endl;
  
  resize(newSize);
  mConsensusQuality.resize(newSize);
  mPaddedBases = size();

  //cout<<"newSize="<<newSize<<endl;
  //cout<<"New C size="<<size()<<endl;
  //cout<<"TempC.size="<<TempC.size()<<endl;

  //cout<<"mConsensusQuality.size="<<mConsensusQuality.size()<<endl;
  //cout<<"TempC.mConsensusQual.size="<<TempC.mConsensusQuality.size()<<endl;

  for(int i = 0;i< newSize;i++){
    (*this)[i] = TempC[i];
    mConsensusQuality[i] = TempC.mConsensusQuality[i];
  }
}

//----------------------------------------
// * TemplateCoverIsZero
//----------------------------------------
// Returns true if the template coverage for 
// this contig is everywhere zero.  False otherwise. 
bool
Contig::TemplateCoverIsZero()
{
  vector<int> TC;
  GetTemplateCoverage(TC);

  bool bIsZero = true;
  for(int i = 0;i< TC.size();i++){
    if (TC[i] > 0) bIsZero = false;
  }

  return(bIsZero);
}


//----------------------------------------
// * MultiSplitContig
//----------------------------------------
// Performs multiple splits of the same contig.  Assumes that 
// split pairs already sorted by leftSplit value. 
// 
void 
AceFile::MultiSplitContig(int MinMiddle,int firstCIdx, vector<int> &leftSplits,vector<int> &rightSplits)
{
  // As I split the contigs, the remaining split points may be in one
  // of the newly created contigs.  SplitContig is used to keep track of 
  // which split goes with which contig.   Initially, of course, everything
  // is in contig cIdx;
  vector<int> SplitContig;
  SplitContig.resize(leftSplits.size());
  for(int i = 0;i< SplitContig.size();i++) SplitContig[i] = firstCIdx;

  // Since deleting contigs effectively renumbers the remaining
  // ones, it is confusing to delete as you go along.  So, instead, 
  // I'm going to keep a list of contigs to delete and do them all 
  // at the end. 
  vector<string> DeleteList;
  
  // Perform each of the splits, unless marked to skip ....
  for(int i = 0; i < leftSplits.size();i++){

    int cIdx = SplitContig[i];	 
    if (cIdx < 0) {
      //cout<<"Skipping split."<<endl;
      continue;  // -1 marks splits to skip. 	 
    }
    Contig &C = (*this)[cIdx]; 	 
    int CSize0 = C.size();
	 
	 
    Contig NewC1,NewC2;
    C.SplitContig(leftSplits[i],rightSplits[i],NewC1,NewC2);
	 
    cout<<"\tC.size="<<C.size()<<"\tNewC1.size="<<NewC1.size()<<"\tNewC2.size()="<<NewC2.size()<<endl;
    cout<<"\tCreads="<<C.mReads.size()<<"\tNewC1reads="<<NewC1.mReads.size()<<"\tNewC2Reads="<<NewC2.mReads.size()<<endl;
    //cout<<"Min Middle="<<MinMiddle<<endl;

    // Default values indicate to AdjustSplitPoints that these contigs are empty.
    int newCIdx=-1,newCIdx2=-1;

    // KJD The reference C is unstable to changes in the vector. 
    // KJD If the resize() commands for NewC1 and NewC2 happen before
    // KJD this reference, C may no longer be valid.  
    // C is already there, so add it to the delete list if it has 
    // zero template coverage and size less than the minimum. 
    if (C.size() > 0){
      if ((C.TemplateCoverIsZero()) && (C.size() < MinMiddle)){
	DeleteList.push_back(C.mName);
      }
    }else{
      // Delete it if it is completely empty. 
      DeleteList.push_back(C.mName);
    }


    // Only add pieces with non-zero template coverage or which are 
    // bigger than minMiddle
    if (NewC1.size() > 0){
      if ((!NewC1.TemplateCoverIsZero()) || (NewC1.size() >= MinMiddle)){
	cout<<"\tAdd C1 "<<NewC1.mName<<endl;
	push_back(NewC1);
	newCIdx = size() -1;
      }  
    }
	 
    if (NewC2.size() > 0){
      if ((!NewC2.TemplateCoverIsZero()) || (NewC2.size() >= MinMiddle)){
	cout<<"\tAdd C2 "<<NewC2.mName<<endl;
	push_back(NewC2);
	newCIdx2 = size() -1;
      }
    }

    CleanInitFromContigs();
	 
    // Adjust Split Points
    AdjustSplitPoints(i,firstCIdx,newCIdx,newCIdx2,leftSplits,rightSplits,SplitContig);	 
  }

  // Now go through and delete the contigs that need to be 
  // deleted. They are being deleted by name to ensure that there
  // are no renumbering problems. 
	
  //	cout<<"BEFORE DELETE size="<<size()<<" DeleteList="<<DeleteList.size()<<endl;

  for(int i = 0;i< DeleteList.size();i++){
    cout<<"\tDeleting "<<DeleteList[i]<<"..."<<flush;
    int deleteIdx = CName2CIdx(DeleteList[i]);
	 
    if ((deleteIdx < 0) || (deleteIdx >= size())){
      cerr<<"ERROR: Trying to delete contig "<<DeleteList[i]<<" "<<deleteIdx<<" out of "<<size()<<endl;
      continue;
    }

    erase(begin()+deleteIdx);
    CleanInitFromContigs(); // KJD This is not necessary. 
    cout<<"done."<<endl;
  }


  //cout<<"AFTER DELETE  size="<<size()<<endl;

}

class ReadSize{
 public:
  int mCIdx;
  int mReadSize;
  
  friend bool operator<(const ReadSize &a,const ReadSize &b){
    return(b.mReadSize < a.mReadSize);
  }
  
};


//----------------------------------------
// * AdjustSplitPoints
//----------------------------------------
// Part of MultiSplitContig.  Once a contig has been split, all of the 
// remaining split points need to be assigned to one of the newly created 
// contigs.  The biggest contig which contains both split points is the 
// one preferred.   If no contig contains both split points, the split 
// points are marked to skip. 
// 
void 
AceFile::AdjustSplitPoints(int currentPair, int firstCIdx,int newCIdx,int newCIdx2,
			   vector<int> &leftSplits,vector<int> &rightSplits,
			   vector<int> &SplitContig)
{
  
  vector<ReadSize> ReadSizes;
  ReadSizes.resize(3);
  ReadSizes[0].mCIdx = firstCIdx;
  ReadSizes[1].mCIdx = newCIdx;
  ReadSizes[2].mCIdx = newCIdx2;

  // Some of the contigs may not exist, so must check for that case. 
  if ((firstCIdx >= 0) && (firstCIdx < size())){
    ReadSizes[0].mReadSize = (*this)[firstCIdx].mReads.size();
  }else {
    ReadSizes[0].mReadSize = 0;
  }

  if ((newCIdx >= 0) && (newCIdx < size())){
    ReadSizes[0].mReadSize = (*this)[newCIdx].mReads.size();
  }else {
    ReadSizes[0].mReadSize = 0;
  }
  
  if ((newCIdx2 >= 0) && (newCIdx2 < size())){
    ReadSizes[0].mReadSize = (*this)[newCIdx2].mReads.size();
  }else {
    ReadSizes[0].mReadSize = 0;
  }

  sort(ReadSizes.begin(),ReadSizes.end());

  // Figure out where each of the new splits goes...
  for(int j = currentPair+1;j< leftSplits.size();j++){		
    int ls = leftSplits[j];
    int rs = rightSplits[j];
	 
    // Find the largest (most reads) contig that contains both ls and rs
    bool bFound = false;
    int  MatchContig = -1;
    for(int k = 0;k < ReadSizes.size(); k++){
      // Skip non-existent contigs. 
      if ((ReadSizes[k].mCIdx >= 0) && (ReadSizes[k].mCIdx < size())){
	Contig &TestC = (*this)[ReadSizes[k].mCIdx];
		
	int ltest = ls - TestC.mOldLeftMost;
	int rtest = ls - TestC.mOldLeftMost;
		
	if ((ltest >=0) && (rtest < TestC.size())){
	  bFound = true;
	  MatchContig = k;
	  break;
	}
      }
    }
	 
    if (bFound){
      int foundCIdx = ReadSizes[MatchContig].mCIdx;
      SplitContig[j] = foundCIdx;
      leftSplits[j] = leftSplits[j] - (*this)[foundCIdx].mOldLeftMost;
      rightSplits[j] = rightSplits[j] - (*this)[foundCIdx].mOldLeftMost;
      //		cout<<"AdjustSplitPoints: Ace.size="<<size()<<"matchContig="<<MatchContig<<" foundCIdx="<<foundCIdx;
      //cout<<"left="<<leftSplits[j]<<"right="<<rightSplits[j]<<endl;
		
    }else{
      SplitContig[j] = -1;		  
    }	
  }
}


//----------------------------------------
// * SplitContig
//----------------------------------------
// Splits the contig at the left and right split points.  Reads
// which overlap these points are retained if they have a mate pair
// further in the contig, and are dropped otherwise. 
//
// Remember: once you split a contig you will need to do CleanInitFromContigs
// on the AceFile to update datastructures there.  
// 
void
Contig::SplitContig(int leftSplit, int rightSplit,Contig &NewC1,Contig &NewC2){
  	
  vector<string> LeftVec,RightVec,MiddleVec;

  FindReadSets(leftSplit, rightSplit, 
	       LeftVec,MiddleVec,RightVec);

  cout<<"\t"<<mName<<"\tSplitContig: "<<leftSplit<<"\t"<<rightSplit;
  cout<<"\t"<<LeftVec.size()<<"\t"<<MiddleVec.size()<<"\t"<<RightVec.size()<<endl;

  int leftPoint,rightPoint;

  // Split off the rightmost reads first and put them into NewC2
  if (RightVec.size() > 0){
    NewC2.mName = mName+"R";
    SplitContigRight(RightVec,NewC2);
  }

  // Then, working leftward, split off the middle reads and put them into NewC1
  if (MiddleVec.size() > 0){
    NewC1.mName = mName+"M";
    SplitContigRight(MiddleVec,NewC1);
  }

  mName=mName+"L";
}

void 
Contig::FindConsensusChangePoints(vector<string> &readList,int &leftmost,int &rightmost)
{
  leftmost = size();
  rightmost = 0;

  for(int i = 0;i< readList.size();i++){
    string rName = readList[i];
    if (mReadIdxMap.Contains(rName)){
      int rIdx = mReadIdxMap[rName];
      int leftIdx = mAlignStart[rIdx];
      int rightIdx = mAlignStart[rIdx] + mReads[rIdx].size() -1;
      if (leftIdx < leftmost) leftmost = leftIdx;
      if (rightIdx > rightmost) rightmost = rightIdx;		
    }else{
      cerr<<"ERROR: mReadIdxMap does not contain "<<rName<<endl;
    }
  }

  // Bound it within the contig. 
  if (leftmost < 0) leftmost = 0;
  if (rightmost < 0) rightmost = 0;
  if (leftmost > size()-1) leftmost = size()-1;
  if (rightmost > size()-1) rightmost = size()-1;


}

//----------------------------------------
// * ReverseComplementContig
//----------------------------------------
// Reverse complement contig and all dependent data. 
// 
// Note that a Contig is a kind of Sequence, and 
// Sequence has a ReverseComplement() method.  Be 
// careful not to confuse that method, which comptes 
// the RC of the sequence itself (the consensus sequence
// here), with this method which generates a contig which 
// is the reverse complement of the original, reads and 
// all. 
//
// In specific:
// 
//  AF lines must be updated (in specific C/U designations)
//  Reads must be reverse complemented. 
//  Reads quality must be reversed. 
//  Consensus must be reverse complemented. 
//  Consensus quality must be reversed. 
//  
//  QA lines must be updated. 
//  BS lines must be updated. 
void 
Contig::ReverseComplementContig()
{
  // Reverse complement the contig sequence and contig quality. 
  ReverseComplement();


  //	cout<<"BEFORE"<<endl;
  //for(int i = 0; i < 100;i++){
  //	cout<<mConsensusQuality[i]<<" ";
  //}
  //cout<<endl;

  mConsensusQuality.Reverse();

  //cout<<"AFTER"<<endl;
  //for(int i = 0; i < 100;i++){
  //	cout<<mConsensusQuality[i]<<" ";
  //}
  //cout<<endl;


  // mTrimLeft
  // mTrimRight
	
  int cSize = size();
  for(int i = 0;i < mReads.size();i++){
    mReads[i].ReverseComplement();

    if (mReadsQuality.size() <= i){
      cerr<<"WARNING: There are "<<mReadsQuality.size()<<" read qualities and "<<mReads.size()<<" reads."<<endl;
    }else{
      if (mReadsQuality[i].size() == 0){
	cerr<<"WARNING: "<<mReads[i].mName<<" has no quality in RCContig()."<<endl;
      }else{
	mReadsQuality[i].Reverse();
      }
    }		

    mAlignStart[i] = cSize-mAlignStart[i]+1-mReads[i].size()+1;
		
    // mReadSense is how the read in the original FA file relates to the 
    // the read in this file.  So whatever sense it was before RC-Contig, 
    // it will be the opposite sense now. 
    if (mReadSense[i] == 'C') mReadSense[i] = 'U';
    else mReadSense[i]='C';


    // mReadTrimLeft
    // mReadTrimRight
    // mBaseSegments				
  }
	
  // Recompute mTInfo

  // Right now, just generate a synthetic BS for the 
  // Could recompute what is already there...
  ComputeSyntheticBS();
 	
}



//----------------------------------------
// * MergeContigRight
//----------------------------------------
// Merges contig C into the right side of this 
// contig, with overlap consensus bases overlapping. 
// This function does NOT do any checking to see that 
// the merge makes sense, it just does it.  
//
// Note: StartLocation is taken to be a zero-based 
// contig coordinate, not the Phrap one based coordinate. 
// 
// Note: This function does NOT delete the original contig. 
// This is left up to the calling program. 
//
void
Contig::MergeContigRight(Contig &C,int StartLocation)
{
  // Add C's reads to this contig...
  for(int i = 0;i < C.mReads.size();i++){
    // Fixup the align-start for the new contig. 
    C.mAlignStart[i] = C.mAlignStart[i]+StartLocation;
    C.CopyRead((*this),C.mReads[i].mName);
  }
	
  // Compute the new right-hand endpoint of the contig
  int newRight = 0;
  for(int r = 0;r< mReads.size();r++){
    int nend = RIdx2CIdx(mReads[r].size()-1,r);
    if (nend > newRight) newRight = nend;
  }

  // Make space for the new consensus. 
  resize(newRight+1);
  mConsensusQuality.resize(newRight+1);

  // Compute a new consensus for overlap position rightward. 
  // KJD - what to do about reads that overlap the StartLocation????
  ComputeConsensus(StartLocation,newRight);

  mPaddedBases = size();
  mNumEntries = mReads.size();
  ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();
	
}

//----------------------------------------
// * AddReads
//----------------------------------------
// Add some reads and associated read qualities to the contig
// at the given contig offsets. Input offsets are zero based. 
// They will be converted to one-based when writing AceFile.
// 
void
Contig::AddReads(SequenceVector &Reads,SequenceQualityVector &ReadsQual,vector<int> &offsets)
{
  int oldSize = mReads.size();	
  int newSize = oldSize+Reads.size();
  int nextOffset = oldSize;


  mReads.resize(newSize);
	
  mReadsQuality.resize(newSize);

  mAlignStart.resize(newSize);
  mReadSense.resize(newSize);

  for(int i = 0;i< Reads.size();i++){
    mReads[nextOffset] = Reads[i];
    mReadsQuality[nextOffset] = ReadsQual[i];

    // Offsets are presumed to be
    mAlignStart[nextOffset] = offsets[i]+1;
    mReadSense[nextOffset] = 'U';
    nextOffset++;
  }

  mNumEntries = mReads.size();
  ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();

  //  vector<int>           mReadTrimLeft;
  //vector<int>           mReadTrimRight;
 
}





//----------------------------------------
// * AddReads
//----------------------------------------
// Add some reads to the contig at the given contig offsets. 
// 
// KJD HORRIBLE HORRIBLE  The other version of AddReads adds one 
// KJD HORRIBLE HORRIBLE  to each offset value because the calling 
// KJD                    function doesn't.  However, the version that 
// KJD                    calls this version does add one already.  
// KJD                    This MUST change as soon as the crunch is over. 
// 
void
Contig::AddReads(SequenceVector &Reads,vector<int> &offsets)
{
  int oldSize = mReads.size();	
  int newSize = oldSize+Reads.size();
  int nextOffset = oldSize;
  mReads.resize(newSize);
  mAlignStart.resize(newSize);
  mReadSense.resize(newSize);
	
  if ((offsets.size() != Reads.size()) || (offsets.size() ==0)) {
    cerr<<"AddReads ERROR: offset.size="<<offsets.size()<<" Reads.size="<<Reads.size()<<endl;
    return;
  }

  for(int i = 0;i< Reads.size();i++){
    mReads[nextOffset] = Reads[i];
    mAlignStart[nextOffset] = offsets[i];
    mReadSense[nextOffset] = 'U';
    nextOffset++;
  }

  mNumEntries = mReads.size();
  ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();

  //  vector<int>           mReadTrimLeft;
  //vector<int>           mReadTrimRight;
 
}



//----------------------------------------
// * AddReads
//----------------------------------------
// Add some reads to the contig at the given contig offsets. 
// 
// KJD HORRIBLE HORRIBLE  The other version of AddReads adds one 
// KJD HORRIBLE HORRIBLE  to each offset value because the calling 
// KJD                    function doesn't.  However, the version that 
// KJD                    calls this version does add one already.  
// KJD                    This MUST change as soon as the crunch is over. 
// 
void
Contig::AddReadsNoBS(SequenceVector &Reads,vector<int> &offsets,vector<char> &Senses)
{
  int oldSize = mReads.size();	
  int newSize = oldSize+Reads.size();
  int nextOffset = oldSize;
  mReads.resize(newSize);
  mAlignStart.resize(newSize);
  mReadSense.resize(newSize);
	
  if ((offsets.size() != Reads.size()) || (offsets.size() ==0)) {
    cerr<<"AddReads ERROR: offset.size="<<offsets.size()<<" Reads.size="<<Reads.size()<<endl;
    return;
  }

  for(int i = 0;i< Reads.size();i++){
    mReads[nextOffset] = Reads[i];
    mAlignStart[nextOffset] = offsets[i];
    mReadSense[nextOffset] = Senses[i];
    nextOffset++;
  }

  mNumEntries = mReads.size();
  //	ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();

  //  vector<int>           mReadTrimLeft;
  //vector<int>           mReadTrimRight;
}



//----------------------------------------
// * AddReads
//----------------------------------------
// Add some reads to the contig at the given contig offsets. 
// 
// KJD HORRIBLE HORRIBLE  The other version of AddReads adds one 
// KJD HORRIBLE HORRIBLE  to each offset value because the calling 
// KJD                    function doesn't.  However, the version that 
// KJD                    calls this version does add one already.  
// KJD                    This MUST change as soon as the crunch is over. 
// 
void
Contig::AddReadsNoBS(SequenceVector &Reads,vector<int> &offsets)
{
  int oldSize = mReads.size();	
  int newSize = oldSize+Reads.size();
  int nextOffset = oldSize;
  mReads.resize(newSize);
  mAlignStart.resize(newSize);
  mReadSense.resize(newSize);
	
  if ((offsets.size() != Reads.size()) || (offsets.size() ==0)) {
    cerr<<"AddReads ERROR: offset.size="<<offsets.size()<<" Reads.size="<<Reads.size()<<endl;
    return;
  }

  for(int i = 0;i< Reads.size();i++){
    mReads[nextOffset] = Reads[i];
    mAlignStart[nextOffset] = offsets[i];
    mReadSense[nextOffset] = 'U';
    nextOffset++;
  }

  mNumEntries = mReads.size();
  //	ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();

  //  vector<int>           mReadTrimLeft;
  //vector<int>           mReadTrimRight;
 
}


//----------------------------------------
// * SplitContigRight
//----------------------------------------
// Same as 0ther Contig::SplitContigRight except there 
// is no kill set. 
// 
void 
Contig::SplitContigRight(vector<string> &RightVec,Contig &NewC)								  
{
  int leftmostAlteredConsensus,rightmostAlteredConsensus;
  FindConsensusChangePoints(RightVec,leftmostAlteredConsensus,rightmostAlteredConsensus);
  
  // Add the RightVec reads to the new right contig
  for(int i = 0;i< RightVec.size();i++){	 
    CopyRead(NewC,RightVec[i]);
  }

  // The align starts for those reads are relative to the 
  // old coordinates.   We need to fix them up. 
  int oldLeft=0,oldRight=0;
  NewC.FixupAlignStarts(0,oldLeft,oldRight);

  int newRight = 0;
  for(int r = 0;r< NewC.mReads.size();r++){
    int nend = NewC.RIdx2CIdx(NewC.mReads[r].size()-1,r);
    if (nend > newRight) newRight = nend;
  }

  // Copy the old consensus and quality up to the rightmost altered point. 
  //NewC.CopyConsensus(0,C,rightmostAlteredConsensus+1,C.size()-1);
  // Compute a new consensus for the potenti
  // Note that rightmostAlteredConsensus is with respect to the original 
  // contig C.  So we have to adjust the end to fit this contig:
  //  
  //NewC.ComputeConsensus(0,rightmostAlteredConsensus-oldLeft);

  // KJD  There is some bug in the CopyConsensus way of creating a new
  // KJD  contig.  Just computing the consensus from scratch seems to work OK. 

  NewC.resize(newRight+1);
  NewC.mConsensusQuality.resize(newRight+1);

  //cout<<"newRight="<<newRight<<endl;
  NewC.ComputeConsensus(0,newRight);

  NewC.mPaddedBases = NewC.size();
  NewC.mNumEntries = NewC.mReads.size();


  cout<<NewC.mName<<"\tNewC size="<<NewC.size()<<" NumReads: "<<NewC.mReads.size()<<endl;


  NewC.ComputeSyntheticBS();
  NewC.mNumBaseSeqs = NewC.mBaseSegments.size();

  string bugName = "kjdtest/"+NewC.mName;
  ofstream bout;
  bout.open(bugName.c_str());
  NewC.WriteAce(bout);
  bout.close();


  // Delete the moved reads from the original contig.  
  for(int i = 0;i< RightVec.size();i++){
    RemoveRead(RightVec[i]);
  }

  mPaddedBases = size();
  mNumEntries = mReads.size();

  cout<<mName<<"\tsize="<<size()<<" NumReads: "<<mReads.size()<<endl;

  mNumBaseSeqs = mBaseSegments.size();
  
  // Most of this consensus should be OK, just recompute the bit where
  // reads may have been removed. 
  if (size() >0){
    ComputeConsensus(leftmostAlteredConsensus,size()-1);
    TrimPadTail();
  }

  ComputeSyntheticBS();

  cout<<"After compute consensus: "<<mName<<"\tsize="<<size()<<" NumReads: "<<mReads.size()<<endl;

}


//----------------------------------------
// * SplitContigRight
//----------------------------------------
// Splits off the right portion of a contig, making a new contig out of it. 
// The region of split is defined by the set of reads in it.   Reads in 
// the KillVec will be omitted from either the old or the new contig.  Reads
// in the RightVec will be removed from the old contig and put into the new one. 
// All remaining reads will stay in the old contig.   The old and new contigs will
// preserve their original consensus and quality up to the point where the 
// moved or removed reads might have an effect.   A new consensus will be computed
// from that point onward.   The new contig will have the same name as the old 
// one with the letter 'A' attached to the end (e.g. Contig68 will split into 
// Contig68 and Contig68A.  
//
// 
void 
Contig::SplitContigRight( vector<string> &KillVec,
			  vector<string> &RightVec,
			  int leftmostAlteredConsensus,
			  int rightmostAlteredConsensus,
			  Contig &NewC)
{

  // Now, split the right half of this contig off and make a new contig
  // with it.
  //Contig NewC;
  NewC.mName = mName+"A";
  
  // Add the RightVec reads to the new right contig
  for(int i = 0;i< RightVec.size();i++){	 
    CopyRead(NewC,RightVec[i]);
  }

  // The align starts for those reads are relative to the 
  // old coordinates.   We need to fix them up. 
  int oldLeft=0,oldRight=0;
  NewC.FixupAlignStarts(0,oldLeft,oldRight);

  int newRight = 0;
  for(int r = 0;r< NewC.mReads.size();r++){
    int nend = NewC.RIdx2CIdx(NewC.mReads[r].size()-1,r);
    if (nend > newRight) newRight = nend;
  }

  // Copy the old consensus and quality upt to the rightmost altered point. 
  //NewC.CopyConsensus(0,C,rightmostAlteredConsensus+1,C.size()-1);
  // Compute a new consensus for the potenti
  // Note that rightmostAlteredConsensus is with respect to the original 
  // contig C.  So we have to adjust the end to fit this contig:
  //  
  //NewC.ComputeConsensus(0,rightmostAlteredConsensus-oldLeft);

  // KJD  There is some bug in the CopyConsensus way of creating a new
  // KJD  contig.  Just computing the consensus from scratch seems to work OK. 

  NewC.resize(newRight+1);
  NewC.mConsensusQuality.resize(newRight+1);

  cout<<"newRight="<<newRight<<endl;
  NewC.ComputeConsensus(0,newRight);

  //cout<<"rightmostAltered="<<rightmostAlteredConsensus<<"\toldLeft="<<oldLeft<<"\t";
  //cout<<"\toldRight = "<<oldRight<<"\tnewRight="<<newRight<<"\tNewC.size="<<NewC.size()<<endl;

  NewC.mPaddedBases = NewC.size();
  NewC.mNumEntries = NewC.mReads.size();

  NewC.ComputeSyntheticBS();
  NewC.mNumBaseSeqs = NewC.mBaseSegments.size();

  // Delete the moved reads from the original contig.  
  for(int i = 0;i< RightVec.size();i++){
    RemoveRead(RightVec[i]);
  }

  // Remove the Kill reads from the middle, if they are still around. 
  for(int i = 0;i< KillVec.size();i++){	 
    RemoveRead(KillVec[i]);
  }

  mPaddedBases = size();
  mNumEntries = mReads.size();
  
  // Most of this consensus should be OK, just recompute the bit where
  // reads may have been removed. 
  ComputeConsensus(leftmostAlteredConsensus,size()-1);

  ComputeSyntheticBS();
  mNumBaseSeqs = mBaseSegments.size();
  TrimPadTail();
}



//----------------------------------------
// * ComputeConsensus
//----------------------------------------
// Recomputes the consensus for the contig index range specified. 
// Currently the consensus is computed by simple majority.  
// The quality is changed to be the sum of the consensus matching
// base qualities, capped at 90 max.    
// 
void
Contig::ComputeConsensus(int sIdx,int eIdx)
{
  if (sIdx < 0) sIdx = 0;
  if (eIdx > size()-1) eIdx = size()-1;


  //cout<<"Compute consensus "<<mName<<"\t"<<sIdx<<"\t"<<eIdx<<endl;
  for(int cIdx = sIdx;cIdx <=eIdx;cIdx++){
    char consensus='N';
    int quality = 0;
    ComputeConsensus(cIdx,consensus,quality);
    //	 cout<<cIdx<<"\t"<<consensus<<"\t"<<quality<<endl;
    (*this)[cIdx] = consensus;
    (*this).mConsensusQuality[cIdx] = quality;
  }

  // Invalidate the old BS lines. 
  mBaseSegments.clear();

}


//----------------------------------------
// * ComputeConsensus2 (debug test version) 
//----------------------------------------
// Computes the consensus and quality at base position cIdx
// Of the reads in this column, the base with the highest 
// total quality is selected to be the consensus, and the 
// consensus quality is the sum of the qualities of the consensus
// matching bases (capped at 90).  Exceptions to this rule
// occur when pads are present.  If the majority of symbols
// are pads, regardless of quality, the consensus is a pad
// (that is, a '*' character). 
//
void 
Contig::ComputeConsensus2(int cIdx,char &consensus,int &quality)
{

  vector<char> bases;
  vector<int>  qualities;
  
  // Figure out what bases are at this consensus index. 
  for(int r = 0;r< mReads.size();r++){
    int readIdx = CIdx2RIdx(cIdx,r);
    if ((readIdx>=0) && (readIdx<mReads[r].size())){
      bases.push_back(mReads[r][readIdx]);
      qualities.push_back(mReadsQuality[r][readIdx]);
    }
  }

  cout<<"bases.size()="<<bases.size()<<endl;

  // Sanity checking.  This also handles the cases where we have 
  // removed all of the reads from a region.   No reads = no sequence. 
  if (bases.size() == 0){
    consensus='*';
    quality = 0;
    return;
  }

  // Count the As, Cs, etc. at this location. 
  int AScore=0,CScore=0,GScore=0,TScore=0,StarCnt=0,totalQuality=0;
  for(int i = 0;i< bases.size();i++){

    if (toupper(bases[i]) == 'A') AScore+= qualities[i];
    if (toupper(bases[i]) == 'C') CScore+= qualities[i];
    if (toupper(bases[i]) == 'G') GScore+= qualities[i];
    if (toupper(bases[i]) == 'T') TScore+= qualities[i];

    if (toupper(bases[i]) == '*'){
      StarCnt++;
      continue;
    }
    totalQuality +=qualities[i];
  }

  cout<<"totalQualit = "<<totalQuality<<endl;

  // If more than 1/2 of the bases are '*', call it as a '*'
  double StarFrac = (double)StarCnt/(double) bases.size();
  if (StarFrac > 0.5){
    consensus = '*';
    quality = 0;
    return;
  }

  // Compute the quality fraction of each base. 
  //double AFrac = (double)AScore/(double)totalQuality;
  //double CFrac = (double)CScore/(double)totalQuality;
  //double GFrac = (double)GScore/(double)totalQuality;
  //double TFrac = (double)TScore/(double)totalQuality;
  
  // Find the maximum quality fraction. This is the consensus
  int maxScore = 0;
  maxScore = max(AScore,CScore);
  maxScore = max(GScore,maxScore);
  maxScore = max(TScore,maxScore);

  cout<<"maxScore="<<maxScore<<endl;

  int matchCount = 0;
  if (AScore == maxScore){	 
    matchCount++;
    if (AScore >=20){
      consensus = 'A';
    }
  }
  if (CScore == maxScore){ 
    matchCount++;
    if (CScore >=20){
      consensus = 'C';
    }
  }
  if (GScore == maxScore){
    matchCount++;
    if (GScore >=20){
      consensus = 'G';
    }
  }
  if (TScore == maxScore){ 
    matchCount++;
    if (TScore >=20){
      consensus = 'T';
    }
  }

  cout<<"matchCount = "<<matchCount<<"\t"<<consensus<<endl;

  // If more than one base has the same quality fraction, 
  // the base call is ambiguous and must be called an N. 
  if (matchCount > 1) {
    consensus = 'N';
    quality = 0;
  }else{
    // Compute the consensus quality
    quality = 0;
    for(int i = 0;i < qualities.size();i++){
      if (toupper(bases[i]) == consensus) {
	quality+=qualities[i]; 
      }
    }
    cout<<"quality="<<quality<<endl;
  }

  // Limit quality to 90 max and 0 min. 
  if (quality>90) quality=90;
  if (quality<0) quality = 0;
 
}


//----------------------------------------
// * ComputeConsensus
//----------------------------------------
// Computes the consensus and quality at base position cIdx
// Of the reads in this column, the base with the highest 
// total quality is selected to be the consensus, and the 
// consensus quality is the sum of the qualities of the consensus
// matching bases (capped at 90).  Exceptions to this rule
// occur when pads are present.  If the majority of symbols
// are pads, regardless of quality, the consensus is a pad
// (that is, a '*' character). 
//
void 
Contig::ComputeConsensus(int cIdx,char &consensus,int &quality)
{

  vector<char> bases;
  vector<int>  qualities;
  
  // Figure out what bases are at this consensus index. 
  for(int r = 0;r< mReads.size();r++){
    int readIdx = CIdx2RIdx(cIdx,r);
    if ((readIdx>=0) && (readIdx<mReads[r].size())){
      bases.push_back(mReads[r][readIdx]);
      qualities.push_back(mReadsQuality[r][readIdx]);
    }
  }

  // Sanity checking.  This also handles the cases where we have 
  // removed all of the reads from a region.   No reads = no sequence. 
  if (bases.size() == 0){
    consensus='*';
    quality = 0;
    return;
  }

  // Count the As, Cs, etc. at this location. 
  int AScore=0,CScore=0,GScore=0,TScore=0,StarCnt=0,totalQuality=0;
  for(int i = 0;i< bases.size();i++){

    if (toupper(bases[i]) == 'A') AScore+= qualities[i];
    if (toupper(bases[i]) == 'C') CScore+= qualities[i];
    if (toupper(bases[i]) == 'G') GScore+= qualities[i];
    if (toupper(bases[i]) == 'T') TScore+= qualities[i];

    if (toupper(bases[i]) == '*'){
      StarCnt++;
      continue;
    }
    totalQuality +=qualities[i];
  }

  // If more than 1/2 of the bases are '*', call it as a '*'
  double StarFrac = (double)StarCnt/(double) bases.size();
  if (StarFrac > 0.5){
    consensus = '*';
    quality = 0;
    return;
  }

  // Compute the quality fraction of each base. 
  //double AFrac = (double)AScore/(double)totalQuality;
  //double CFrac = (double)CScore/(double)totalQuality;
  //double GFrac = (double)GScore/(double)totalQuality;
  //double TFrac = (double)TScore/(double)totalQuality;
  
  // Find the maximum quality fraction. This is the consensus
  int maxScore = 0;
  maxScore = max(AScore,CScore);
  maxScore = max(GScore,maxScore);
  maxScore = max(TScore,maxScore);

  int matchCount = 0;
  if (AScore == maxScore){	 
    matchCount++;
    consensus = 'A';
  }
  if (CScore == maxScore){ 
    matchCount++;
    consensus = 'C';
  }
  if (GScore == maxScore){
    matchCount++;
    consensus = 'G';
  }
  if (TScore == maxScore){ 
    matchCount++;
    consensus = 'T';
  }

  // If more than one base has the same quality fraction, 
  // the base call is ambiguous and must be called an N. 
  if (matchCount > 1) {
    consensus = 'N';
    quality = 0;
  }else{
    // Compute the consensus quality
    quality = 0;
    for(int i = 0;i < qualities.size();i++){
      if (toupper(bases[i]) == consensus) {
	quality+=qualities[i]; 
      }
    }
  }

  // Limit quality to 90 max and 0 min. 
  if (quality>90) quality=90;
  if (quality<0) quality = 0;


  /*   
  // Some sanity checks, if applicable
  if (cIdx >= size()) return;

  int originalQuality = mConsensusQuality[cIdx];
  // See if original was a pad but we called it something else
  if ((*this)[cIdx] == '*'){
  cerr<<cIdx<<"\tOriginal = * \t base="<<consensus<<"\t";
  cerr<<"\tOriginalQuality = "<<originalQuality<<" Quality="<<quality<<endl;
  for(int i = 0;i< bases.size();i++){
  cerr<<"\t"<<bases[i]<<"\t"<<qualities[i]<<endl;
  }
  }

  // See if original was different from our call
  if (toupper((*this)[cIdx]) != toupper(consensus)){
  cerr<<mName<<"\t"<<cIdx<<"\tOriginal = "<<(*this)[cIdx]<<"\t";

  for(int k = 0;k < 10;k++){
  if (cIdx+k < size()){
  cerr<<(*this)[cIdx+k];
  }
  }
  cout<<"\t";


  cerr<<"base="<<consensus<<"\t";
  cerr<<AScore<<"\t"<<CScore<<"\t"<<GScore<<"\t"<<TScore<<"\tMax="<<maxScore<<"\tmatchCnt="<<matchCount<<"\t";
  cerr<<"OriginalQuality = "<<originalQuality<<" Quality="<<quality<<endl;
  for(int i = 0;i< bases.size();i++){
  cerr<<"\t"<<bases[i]<<"\t"<<qualities[i]<<endl;
  }
  }
  */

 
}

//----------------------------------------
// * CreateConsensusForCopy
//----------------------------------------
// Deprecated. 
// 
void 
Contig::CreateConsensusForCopy(){

  // Out with the old. 
  clear();
  mConsensusQuality.clear();

  // All of the alignStart values for this are with respect to the
  // source contig.   We're essentially starting from scratch, so
  // First determine the bounds for the new consensus.  
  int leftmostAlignStart = 99999999;
  for(int r = 0;r < mReads.size();r++){
    if (mAlignStart[r] < leftmostAlignStart){
      leftmostAlignStart = mAlignStart[r];
    }
  }
  
  // OK, leftmostAlignStart is now the zero base position, 
  // so renumber all of the alignStart locations
  int rightmostConsensus = 0;
  for(int r = 0;r< mReads.size();r++){
    mAlignStart[r] = mAlignStart[r] - leftmostAlignStart;
    int rightIdx = mAlignStart[r]+mReads[r].size()-1;
    if (rightIdx > rightmostConsensus) rightmostConsensus = rightIdx;
  }

  int badBaseCount = 0;
  for(int c = 0;c< rightmostConsensus;c++){
    char base;
    int qual;
    ComputeConsensus(c,base,qual);	 
    push_back(base);	 
    mConsensusQuality.push_back(qual);
  }
  mPaddedBases = size();
  mNumEntries = mReads.size();
  mNumBaseSeqs = mBaseSegments.size();
}


//----------------------------------------
// * CopyRead
//----------------------------------------
// Copies a read and all it's ancillary information to a new 
// contig. 
// 
void
Contig::CopyRead(Contig &NewC,string readName)
{
  if (!mReadIdxMap.Contains(readName)){
    cout<<"Can't move "<<readName<<" Doesn't exist."<<endl;
    return;
  }

  int rIdx = mReadIdxMap[readName];

  int newReadIdx = NewC.mReads.size();


  NewC.mReadIdxMap[readName] = newReadIdx;

  NewC.mReads.resize(newReadIdx+1);
  NewC.mReads[newReadIdx] = mReads[rIdx];

  NewC.mReadsQuality.resize(newReadIdx+1);
  NewC.mReadsQuality[newReadIdx] = mReadsQuality[rIdx];

  NewC.mAlignStart.resize(newReadIdx+1);
  NewC.mAlignStart[newReadIdx] = mAlignStart[rIdx];

  NewC.mReadSense.resize(newReadIdx+1);
  NewC.mReadSense[newReadIdx] = mReadSense[rIdx];

  NewC.mQALines.resize(newReadIdx+1);
  NewC.mQALines[newReadIdx] = mQALines[rIdx];

  NewC.mReadTrimLeft.resize(newReadIdx+1);
  NewC.mReadTrimLeft[newReadIdx] = mReadTrimLeft[rIdx];

  NewC.mReadTrimRight.resize(newReadIdx+1);
  NewC.mReadTrimRight[newReadIdx] = mReadTrimRight[rIdx];
  
  // Add in one dummy base segment, just to not break some parsers
  NewC.mBaseSegments.resize(1);
  NewC.mBaseSegments[0].mStart=0;
  NewC.mBaseSegments[0].mEnd=0;
  NewC.mBaseSegments[0].mName = readName;

  string Template = mRead2TemplateMap[readName];
  NewC.mRead2TemplateMap[readName] = Template;
  NewC.mTemplate2ReadMap.Add(Template,newReadIdx);

  NewC.mTrimLeft = 0;
  NewC.mTrimRight= 0;  // KJD ??

  // Some things I'm skipping for now. 
  //  Int2IntMap   mReadTemplateDir; 
  
  // Now that everything has been copied to the new contig, remove
  // the read from the source contig. 
  //C.RemoveRead(readName);
}


//----------------------------------------
// * GetNumReadsIn
//----------------------------------------
// Returns a vertical stripe of read bases at contig position cIdx
// winStart and winEnd are used to determine how many reads to pad 
// with spaces.  That is, there may be reads in the window that do 
// not appear at cIdx, but vertical padding must reflect the whole 
// window's worth of characters.  
// 
// ........ACTGAGATTACTAGGATAC..................
// ..ACATTAACTGAGATTACTGGGATAGGGCAG.............
// ................................AAGGTACTGGGGG
int 
Contig::GetNumReadsInWindow(int winStart,int winEnd)
{
  string outString;
  int rSize = mReads.size();

  // Examine each read to see if it is in our window and count it.  
  int numReadsInWin = 0;
  for(int r = 0;r< rSize;r++){
    for(int c = winStart;c<winEnd;c++){				
      int rIdx = CIdx2RIdx(c,r);
      if ((rIdx > 0) && (rIdx < mReads[r].size())){
	numReadsInWin++;
	break;
      }
    }
  }
  return(numReadsInWin);
}


//----------------------------------------
// * GetWindowReadName
//----------------------------------------
// Returns the rIdx'th read in the defined window.  
string
Contig::GetWindowReadName(int winStart,int winEnd,int targetIdx)
{
  string outString;
  int rSize = mReads.size();

  IntSet WinReads;

  // Examine each read to see if it is in our window and if so 
  // if it's the rth such read.  
  int readCount = 0;  
  for(int r = 0;r< rSize;r++){
    for(int c = winStart;c<winEnd;c++){				
      int rIdx = CIdx2RIdx(c,r);
      if ((rIdx > 0) && (rIdx < mReads[r].size())){
	if (readCount == targetIdx) return(mReads[r].mName);
	readCount++; // Found another read in our window...
	break;
      }
    }
  }
  return("NOT FOUND");
}



//----------------------------------------
// * GetAlignedCharacters
//----------------------------------------
// Returns a vertical stripe of read bases at contig position cIdx
// winStart and winEnd are used to determine how many reads to pad 
// with spaces.  That is, there may be reads in the window that do 
// not appear at cIdx, but vertical padding must reflect the whole 
// window's worth of characters.  
// 
// ........ACTGAGATTACTAGGATAC..................
// ..ACATTAACTGAGATTACTGGGATAGGGCAG.............
// ................................AAGGTACTGGGGG
string
Contig::GetAlignedCharacters(int winStart,int winEnd,int cIdx)
{
  string outString;
  int rSize = mReads.size();

  // Examine each read to see if it is in our window and if so 
  // if it is also in our column.  If in our window we add a pad
  // character to our output.  If in column we add the read character
  // to out output. 
  int numReadsInWin = 0;
  for(int r = 0;r< rSize;r++){
    for(int c = winStart;c<winEnd;c++){				
      int rIdx = CIdx2RIdx(c,r);
      if ((rIdx > 0) && (rIdx < mReads[r].size())){

	// Read is in the window, see if it is also in the column
	int columnIdx = CIdx2RIdx(cIdx,r);
	if ((columnIdx > 0) && (columnIdx < mReads[r].size())){
	  outString+=mReads[r][columnIdx];
	}else{
	  outString+=" ";
	}
	// Either way, we're done with this read. 
	break;
      }
    }
  }
  return(outString);
}


//----------------------------------------
// * GetAlignedCharacter
//----------------------------------------
//
char
Contig::GetAlignedCharacter(int readNum,int ConsensusPosition)
{
  if (readNum >=mReads.size()) return('#');
  if (ConsensusPosition >= size()) return('@');

  int readIdx = CIdx2RIdx(ConsensusPosition,readNum);

  if ((readIdx >= mReadTrimLeft[readNum]) && (readIdx <= mReadTrimRight[readNum])){
    return(mReads[readNum][readIdx]);
  }else if (readIdx >= mReads[readNum].size()){
    return('.');
  }else{
    if (toupper(mReads[readNum][readIdx] == 'X')){
      return('%'); // Symbol for a trimmed X
    }else{
      return('-');
    }
  }
}


//----------------------------------------
// * Printalignment
//----------------------------------------
//
void
Contig::PrintAlignment(ostream &fout,int start,int end)
{
  Contig &C = (*this);

  if (end < 0) end = C.size() -1;

  int c = 0;
  int lineStart = start;
  do{
    // Print consensus
    //fout.width(20);
    //fout<<C.mName;
    PrintSpacePaddedString(fout,16,C.mName);
    fout.width(7);
    fout<<lineStart<<" ";
    for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
      if ((idx >= C.mTrimLeft) && (idx <= C.mTrimRight)){
	fout<<C[idx];
      }else if (idx > C.size()){
	fout<<".";
      }else{
	fout<<"-";  // quality trimmed. 
      }
    }
    fout<<endl;

    // Print each of the reads. 
    for(int r = 0;r< C.mReads.size();r++){
      int readStart = CIdx2RIdx(0,r);
      int readLineStart = readStart + lineStart;

      PrintSpacePaddedString(fout,16,C.mReads[r].mName);
      fout.width(7);
      fout<<readLineStart<<" ";
      for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
	int readIdx = readStart + idx;
	if ((readIdx >= C.mReadTrimLeft[r]) && (readIdx <= C.mReadTrimRight[r])){
	  fout<<C.mReads[r][readIdx];
	}else if (readIdx >= C.mReads[r].size()){
	  fout<<".";
	}else{
	  if (toupper(C.mReads[r][readIdx] == 'X')){
	    fout<<"%"; // Symbol for a trimmed X
	  }else{
	    fout<<"-"; // Symbol for a trimmed whatever else. 
	  }
	}
      }
      fout<<endl;
    }
  	 
    fout<<endl;
    lineStart+=mPrintLineSize;
  }while((lineStart < C.size()) && (lineStart <=end));
}

//----------------------------------------
// * PrintalignmentWithQuality
//----------------------------------------
//
void
Contig::PrintAlignmentWithQuality(ostream &fout,int start,int end)
{
  Contig &C = (*this);

  if (end < 0) end = C.size() -1;

  int c = 0;
  int lineStart = start;
  do{
    // Print consensus
    //	 fout.width(20);
    //fout<<C.mName;
    PrintSpacePaddedString(fout,16,C.mName);
    fout.width(7);
    fout<<lineStart<<" ";
    for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
      if ((idx >= C.mTrimLeft) && (idx <= C.mTrimRight)){
	fout<<" "<<C[idx]<<" ";
      }else if (idx > C.size()){
	fout<<" . ";
      }else{
	fout<<" - ";  // quality trimmed. 
      }
    }
    fout<<endl;

    // Print consensus quality
    PrintSpacePaddedString(fout,16," ");
    fout.width(7);
    fout<<lineStart<<" ";
    for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
      if ((idx >= C.mTrimLeft) && (idx <= C.mTrimRight)){
	PrintSpacePaddedString(fout,3,itoa(C.mConsensusQuality[idx]));
	//fout<<C.mConsensusQuality[idx]<<"  ";
      }else if (idx > C.size()){
	fout<<".  ";
      }else{
	fout<<"-  ";  // quality trimmed. 
      }
    }
    fout<<endl<<endl;

    // Print each of the reads. 
    for(int r = 0;r< C.mReads.size();r++){
      // Print read sequence. 
      int readStart = -C.mAlignStart[r]+1;   // KJD CHECK mAlignStart is one based. 
      int readLineStart = readStart+lineStart;
      //fout.width(20);
      //fout<<C.mReads[r].mName;
      PrintSpacePaddedString(fout,16,C.mReads[r].mName);
      fout.width(7);
      fout<<readLineStart<<" ";
      for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
	int readIdx = readStart + idx;
	if ((readIdx >= C.mReadTrimLeft[r]) && (readIdx <= C.mReadTrimRight[r])){
	  fout<<" "<<C.mReads[r][readIdx]<<" ";
	}else if (readIdx >= C.mReads[r].size()){
	  fout<<" . ";
	}else{
	  if (toupper(C.mReads[r][readIdx] == 'X')){
	    fout<<" % "; // Symbol for a trimmed X
	  }else{
	    fout<<" - "; // Symbol for a trimmed whatever else. 
	  }
	}
      }
      fout<<endl;
		
      // Print read quality
      readStart = -C.mAlignStart[r]+1;
      readLineStart = readStart+lineStart;
      PrintSpacePaddedString(fout,16+7+1," ");
      //		fout.width(7+16);
      //fout<<readLineStart<<" ";
      for(int idx = lineStart;idx< lineStart+mPrintLineSize;idx++){
	int readIdx = readStart + idx;
	if ((readIdx < 0) || (readIdx >= C.mReadsQuality[r].size())){
	  fout<<"   ";
	}else{
	  if (C.mReadsQuality[r][readIdx] < 10){
	    fout<<" "<<C.mReadsQuality[r][readIdx]<<" ";
	  }else{
	    fout<<C.mReadsQuality[r][readIdx]<<" ";
	  }
	}
      }
      fout<<endl;
    }
	 
    fout<<endl;
    fout<<endl;
    lineStart+=mPrintLineSize;
  }while((lineStart < C.size()) && (lineStart <=end));
}



//=====================

/* 
//int lReadIdx = C.CIdx2RIdx(leftCrash,r);
//int rReadIdx = C.CIdx2RIdx(rightCrash,r);




string RName = C.mReads[r].mName;
int rSize = C.mReads[r].size();

int rEnd = C.RIdx2CIdx(rSize,r);
int rStart = C.RIdx2CIdx(

// if (rReadIdx <= rSize){
//	cout<<"LTSIze rReadIdx="<<rReadIdx<<" <= size="<<rSize<<endl;
//	rLTSize++;
//}


//if (rReadIdx <=0){ 
//	cout<<"LTZero rReadIdx="<<rReadIdx<<" <=  Zero="<<0<<endl;
//	rLTZero++;
//}

bool bHasGoodMate = C.HasGoodMateInContig(RName);
bool bFirmlyLeft = (lReadIdx >= rEnd);
bool bOverlapsLeftCrash = ((0 < lReadIdx) && (lReadIdx < rEnd));
bool bUnanchoredLeft = (bOverlapsLeftCrash && !bHasGoodMate);
	 

bool bFirmlyRight = (rReadIdx <= rStart);
bool bOverlapsRightCrash = ((0 < rReadIdx) && (rReadIdx < rEnd));
bool bUnanchoredRight = (bOverlapsRightCrash && (!bHasGoodMate));
bool bAnchoredRight = (bOverlapsRightCrash && bHasGoodMate);

bool bBetweenLeftAndRightCrash = ((lReadIdx < 0) && (rReadIdx >= rEnd));

cout<<lReadIdx<<"\t"<<rSize<<"\t";
cout<<"bGoodMate="<<bHasGoodMate<<"\t";
cout<<"bFLeft="<<bFirmlyLeft<<"\t";
cout<<"bOverlapsLeft="<<bOverlapsLeftCrash<<"\t";
cout<<"bUnanchoredLeft="<<bUnanchoredLeft<<endl;

cout<<rReadIdx<<"\t"<<rSize<<"\t";
cout<<"bGoodMate="<<bHasGoodMate<<"\t";
cout<<"bFRight="<<bFirmlyRight<<"\t";
cout<<"bOverlapsRight="<<bOverlapsRightCrash<<"\t";
cout<<"bUnanchoredRight="<<bUnanchoredRight<<endl;
cout<<endl;
	 
if (bFirmlyRight){
MoveSet.push_back(RName);
}else{
if (bBetweenLeftAndRightCrash){
KillSet.push_back(RName);
}
}
*/



//------------------------------------------------------------------------------
// Revision History
//==============================================================================
// 9/24/2002 
// GetReadsQuality read with concatenation, and didn't notice if the same read name
// had already been read.   If the read was duplicated, it read a quality that was 2x 
// the size.  Changed to read quality only if previous quality size is zero.  
// 
// 9/25/2002  
//  Contig::Copy wasn't copying mReadTrimRight and mReadTrimLeft appropriately. 
//      


